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ABSTRACT 

We analyse the non-linear propagation and dissipation of axisymmetric waves in ac- 
cretion discs using the ZEUS-2D hydrodynamics code. The waves are numerically 
resolved in the vertical and radial directions. Both vertically isothermal and thermally 
stratified accretion discs are considered. The waves are generated by means of reso- 
nant forcing and several forms of forcing are considered. Compressional motions are 
taken to be locally adiabatic (7 = 5/3). Prior to non-linear dissipation, the numerical 
results are in excellent agreement with the linear theory of wave channelling in pre- 
dicting the types of modes that are excited, the energy flux by carried by each mode, 
and the vertical wave energy distribution as a function of radius. In all cases, waves 
are excited that propagate on both sides of the resonance (inwards and outwards). 
For vertically isothermal discs, non- linear dissipation occurs primarily through shocks 
that result from the classical steepening of acoustic waves. For discs that are sub- 
stantially thermally stratified, wave channelling is the primary mechanism for shock 
generation. Wave channelling boosts the Mach number of the wave by vertically con- 
fining the wave to a small cool region at the base of the disc atmosphere. In general, 
outwardly propagating waves with Mach numbers near resonance A4 V >0.0T undergo 
shocks within a distance of order the resonance radius. 
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1 INTRODUCTION 

Gaseous discs play a significant role in the evolution of bi- 
nary star and planetary systems. Generally, the tidal forces 
from stellar or planetary companions distort the discs from 
an axisymmetric form. However, where resonances occur in 
the discs, the tidal forces also generate waves that trans- 
port energy and angular momentum. The resulting reso- 
nant torques cause orbital evolution of the perturbing ob- 
jects (e.g. Goldreich & Tremaine 1980; Lin & Papaloizou 
1993; Lubow & Artymowicz 1996). The waves also cause the 
discs to evolve since they transfer their energy and angular 
momentum to the disc in the locations where they damp. 
Damping may be provided by shocks, radiative damping 
(Cassen & Woolum 1996), or other viscous damping mech- 
anisms. 

Waves in discs have also been considered as possible ex- 
planations of quasi-periodic variability in the X-ray emission 
of systems involving accretion discs around black holes (e.g. 
Kato & Fukue 1980; Nowak & Wagoner 1991; Kato 2001). 
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Initial studies of wave propagation approximated the 
disc as two-dimensional and ignored the effects of vertical 
structure (perpendicular to the disc plane). Goldreich & 
Tremaine (1979) developed a two-dimensional linear theory 
for resonant tidal torques and the associated wave propa- 
gation. Later studies of the non-linear case found that the 
torques were within a few percent of those predicted by 
the two-dimensional linear formula (Shu, Yuan & Lissauer 
1985; Yuan & Cassen 1994). For a thin disc that is verti- 
cally isothermal and has an isothermal thermodynamic re- 
sponse (7 = 1), the only wave excited in the disc has a two- 
dimensional structure. The wavefronts are perpendicular to 
the disc plane and the motion is purely horizontal (paral- 
lel to the disc plane). Thus, a two-dimensional treatment 
is valid in a thin disc if both the vertical structure of the 
disc and its thermodynamic response are locally isothermal. 
However, this is not realistic for most discs, such as those in 
cataclysmic variables, circumstellar and circumbinary discs 
around pre-main-sequence stars, and protoplanetary discs. 
In such discs, the waves are no longer two-dimensional and 
the vertical structure of the disc must be taken into account. 

A study of three-dimensional wave propagation was 
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made using linearised numerical simulations (Lin, Pa- 
paloizou & Savonije 1990a, b). Unfortunately, the limited 
range of physical parameter space that was covered led to 
the interpretation that waves in a vertically thermally strat- 
ified disc are refracted upwards into the atmosphere where 
they dissipate via shocks. Later, it was realised that the 
linear problem could be solved semi-analytically (Lubow & 
Pringle 1993; Korycansky & Pringle 1995; Lubow & Ogilvie 
1998; Ogilvie & Lubow 1999). These studies found that the 
propagation of waves in a thin disc is similar to the propa- 
gation of electromagnetic waves along a waveguide (Jackson 
1962; Landau & Lifshitz 1960). The motions in the disc can 
be described in terms of modes. For each mode, the vertical 
wave structure is determined in detail at each disc radius. 
The horizontal variations of the mode are treated by means 
of a WKB approximation in which the horizontal (radial) 
wavenumber is determined as a function of wave frequency. 
A flux conservation law determines the mode amplitude as 
a function of radius. Depending on the vertical structure of 
the disc and the thermodynamic response of the gas, the 
wave energy may be channelled towards the surface of the 
disc, towards the mid-plane of the disc, or occupy the entire 
vertical extent as it propagates in radius. 

Ideally, one would like to determine the non-linear, non- 
axisymmetric response of a variety of disc models (verti- 
cally isothermal and thermally stratified) to resonant forc- 
ing. Unfortunately, this is a formidable numerical problem 
at present. Instead we begin, in this paper, with an anal- 
ysis of the non-linear axisymmetric response. The axisym- 
metric case offers the major advantage of being expressed 
as a two-dimensional problem, which is currently accessible 
numerically. Although axisymmetric waves only transport 
energy and not angular momentum, they should resemble 
the response of the disc to low-azimuthal number tidal forc- 
ing, as arises in binary star systems. The local physics of 
the waveguide is determined within a region whose size is 
a few times the semi-thickness of the disc. On this scale a 
non-axisymmetric wave of low azimuthal wavenumber is al- 
most indistinguishable from an axisymmetric wave. This ex- 
plains why the local physics of the waveguide can be studied 
within the shearing-sheet approximation (Lubow & Pringle 
1993) and why the azimuthal wavenumber appears in the 
dispersion relation only through a Doppler shift of the wave 
frequency. 

Consequently, in this paper, we study the non-linear 
propagation and dissipation of axisymmetric, but fully re- 
solved, waves in accretion discs. The purpose of the paper is 
threefold. First, we wish to determine, using non-linear hy- 
drodynamic calculations, whether or not the semi-analytical 
linear theory provides an accurate description of wave prop- 
agation in accretion discs. Secondly, although linear theory 
can predict where non-linearity is likely to become impor- 
tant, shocks and the process of energy deposition (and an- 
gular momentum deposition with non-axisymmetric waves) 
cannot be handled. We aim to determine how and where 
dissipation occurs through shocks and how accurately this 
can be predicted from the linear analysis. Finally, we wish 
to determine how well non-linear hydrodynamic calculations 
can model the problem and what resolution is required. 

The outline of this paper is as follows. In Section 2, we 
briefly review the ZEUS-2D hydrodynamic code that is used 
to obtain the non-linear results. In Section 3, we discuss the 



structure of the equilibrium discs that we model and the 
grid layout for the ZEUS-2D calculations. Section 4 briefly 
reviews the semi-analytical linear method for solving the 
three-dimensional wave propagation problem, describes our 
method of wave excitation, and discusses the requirements 
for convergence of the numerical calculations. Sections 5 and 
6 contain the results for a vertically isothermal disc, and a 
polytropic disc with a vertically isothermal atmosphere, re- 
spectively. Intermediate discs with optical depths not much 
larger than unity are considered in Section 7. A summary of 
the results is contained in Section 8. 



2 NUMERICAL METHOD 

The non-linear hydrodynamic calculations presented here 
were performed with the ZEUS-2D code developed by Stone 
& Norman (1992). No modifications to the code were re- 
quired, except the addition of damping boundary conditions 
for test purposes only (see below). 

The ZEUS-2D code comes with various options for the 
interpolation scheme and artificial viscosity. We use the 
van Leer (second-order) interpolation scheme. Linear and 
quadratic artificial viscosities are provided in the code, with 
two forms of the quadratic viscosity: the standard von Neu- 
mann & Richtmyer (1950) form, and a tensor form. We use 
only a quadratic artificial viscosity to handle shocks. Calcu- 
lations were performed both with the von Neumann & Richt- 
myer form and with the tensor form. We find no significant 
difference in the results. All the results presented here use 
the standard von Neumann & Richtmyer form. The strength 
of the quadratic viscosity is controlled by a parameter g v isc- 
Generally, we use g v i sc = 4, but we present some results from 
calculations that were performed with g v i sc = 0.5. No other 
type of viscosity (e.g. one that includes a shear stress) is 
used. 

All of the calculations presented here use reflecting 
boundaries. In order to verify that our results were not af- 
fected by reflections from the boundaries, some calculations 
were performed with damping boundaries at the inner and 
outer radii of the numerical grid. Damping is provided by 
adding an acceleration 



to the inner or outer 10 — 20 radial zones of the grid, where v 
is the velocity of the gas and ui is the frequency of the forcing 
that is used to generate the waves (see Section ^) . There is 
no significant difference between the results obtained using 
reflecting and damping boundaries. 



3 MODELLING THE ACCRETION DISCS 

We consider the propagation and dissipation of axisymmet- 
ric (m = 0) waves in two types of accretion disc: a locally 
vertically isothermal disc, and a polytropic disc (polytropic 
index n — 1.5) that joins smoothly on to an isothermal at- 
mosphere. We ignore all effects of disc turbulence. A full de- 
scription of the equilibrium structures is given in Appendix 
A; we only briefly list their properties here. The disc is de- 
scribed in cylindrical coordinates (R,<f),z). As described in 
Appendix A, we shall work in dimensional units, such that 
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the disc angular velocity is unity when the dimensionless 
disc radius R = 1. The discs are Keplerian, with angular 
velocity Q = R~ 3 ^ 2 , apart from small corrections. 

In both discs, the value of the polytropic exponent 
that describes the adiabatic pressure-density relation for the 
waves is 7 = 5/3. For convenience in modelling, we choose 
discs with finite inner and outer radii, Ri and 7?2- We also 
choose a scale-height that increases linearly with radius, 
which requires that the sound speed varies as c oc _R _1//2 . 
The scale-height of the vertically isothermal disc is H oc R. 
The polytropic disc has a well-defined polytropic layer with 
a semi-thickness H p , that smoothly changes into an isother- 
mal atmosphere with its own scale-height Hi, above z w H p . 
Both Hp and Hi are proportional to R. From this point on, 
we will refer to H p simply as H for the polytropic disc. We 
use discs with H/R = 0.05, 0.1, and 0.2. 

The polytropic disc is a fair representation of an opti- 
cally thick disc in which the temperature is greatest at the 
mid-plane and declines with increasing height until the pho- 
tosphere is reached. Above this height, the atmosphere is 
nearly isothermal. We quote an approximate effective opti- 
cal depth at the mid-plane, r, through the relation (e.g. Bell 
et al. 1997) 



(2) 



where T m and Tj are the temperatures at the mid-plane and 
in the isothermal atmosphere, respectively. For our standard 
parameters (Appendix A and Section^), the effective optical 
depth is t « 25000. 

The mid-plane density profiles are power laws, 
p(R,(/>,z — 0) oc i?~ 3 / 2 , throughout the main body of the 
disc. The density drops smoothly at the edges of the discs 
over a distance equal to the local value of H . For the ver- 
tically isothermal discs, we adopt Ri = 0.5 and Ri = 10.0. 
The polytropic discs have outer radii of R2 = 3.0. The value 
of the inner radius is Ri = 0.2 for all but the highest- 
resolution calculations, for which Ri = 0.6. 

To model axisymmetric waves in these discs in ZEUS- 
2D, we use spherical polar coordinates (r,8). For the verti- 
cally isothermal discs, we adopt inner and outer radial grid 
boundaries n n = 0.35, r out = 13, and the grid encompasses 
8 vertical scale-heights above and below the mid-plane (i.e. 
6» min = tt/2 - 8H/r and <9 max = tt/2 + 8H/r). 

For the polytropic discs, the outer grid radius is chosen 
as r out = 4, and 3 vertical scale-heights are modelled above 
and below the mid-plane (except in Section fr| where 6 ver- 
tical scale-heights are modelled). The inner grid radius r ln 
depends on the resolution. We adopt n n = 0.15 for all but 
the highest resolution, for which r- m = 0.6. 

For all cases, we adopt a logarithmic radial grid in which 



(Ar) i+1 /(Ar) ; = V°ut/nn, 



(3) 



where N r is the number of radial grid zones. The 8 grid is 
uniform with 



Ng 



'max i/niin 



0.1 



(H/r) (Ar)j+i/(Ar)i - 1 



(4) 



zones. This choice means that discs with H/r = 0.1 have 
grid zones that have nearly equal radial and vertical lengths. 
Discs with other values of H/r have the same number of 



zones per vertical scale-height and the same radial spacing, 
but each zone is either stretched or compressed vertically. 

In practice, calculations of vertically isothermal discs 
were performed with N r x Ng zones of: 250 x 111, 500 x 221, 
or 1000 x 441. Calculations of the standard polytropic discs 
were computed with: 250 x 45, 500 x 91, 1000 x 183, or 
1156 x 365 zones (the last of which has a resolution equiva- 
lent to 2000 x 365 zones but models a smaller range of radii 
than the lower-resolution calculations). The highest resolu- 
tion calculations took up to 6000 CPU hours on 440 MHz 
Sun Workstations with each factor of 2 increase in resolu- 
tion taking a factor of « 8 longer (a factor of 4 from the 
increase in the number of zones and a factor of 2 from the 
decrease in the time step). In some cases even higher reso- 
lution would be desirable, but doubling the resolution again 
is impractical. 

The models are run until any transient structure disap- 
pears. Typically, for discs with H/r — 0.1, this takes « 30 
orbital periods at the resonant radius (usually r = 1, see 
below). The number of periods required is inversely propor- 
tional to H/r. We calculate the wave energy, mean Mach 
number, radial energy flux, and dissipation rate by averag- 
ing these quantities over N ~ 10 wave periods, P, after the 
disc has completed 40 rotations at the resonant radius. Let 
(v r ,vg) denote the velocity deviations from Keplerian. The 
mean (time-averaged) local wave kinetic energy density is 
calculated as 



E{r,6) 



1 

NR 



-p(v r + Vg) dt. 



(5) 



The mean local Mach number is calculated as 
1 f NP 

M(r,e) = —J ^/ P (v?, + vl)/{ 1P ) dt. (6) 

The mean vertically integrated radial energy flux is calcu- 
lated as 



/r(r) = 27T 



1 

NR 



v r (p — po) dt r sm8 d8 (7) 



where po is the mean pressure, which is calculated over the 
previous wave period. The dissipation rate is calculated as 
the average over several wave periods of the rate, per unit 
volume, at which thermal energy is deposited by the artificial 
viscosity. 



4 AXISYMMETRIC WAVE EXCITATION AND 
PROPAGATION 

The radial propagation of linear axisymmetric waves in a 
vertically isothermal accretion disc was analysed by Lubow 
& Pringle (1993). Subsequently, this work was extended to 
study the radial propagation of waves in vertically ther- 
mally stratified discs (Korycansky & Pringle 1995; Lubow & 
Ogilvie 1998; Ogilvie & Lubow 1999) and magnetized discs 
(Ogilvie 1998). 

4.1 The two-dimensional wave 

The simplest wave structure is that of a two-dimensional 
wave that propagates radially through a vertically isother- 
mal disc and has no vertical motion (e.g. Goldreich & 
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Tremaine 1979). The dispersion relation for the axisymmet- 
ric two-dimensional wave is 



re 2 + ^fc 2 
P 



(8) 



where uj > is the wave frequency, re is the local epicyclic 
frequency of the disc, and k is the radial wavenumber. For 
a Keplerian disc, we have that re = fi and the radius where 
uj = re is an outer Lindblad resonance. Its behaviour in 
launching waves is essentially like that of an outer Lindblad 
resonance encountered in non-axisymmetric cases (e.g. Gol- 
dreich & Tremaine 1979). The two-dimensional axisymmet- 
ric wave can propagate only outside this radius (i.e. where 
uj > fl) for a Keplerian disc. 

If such a wave behaves purely isothermally (i.e. 7 = 1), 
it will propagate without loss to the outer radius of the disc. 
However, as we will see in this paper, adiabatic waves with 
7 = 5/3 steepen into shocks after propagating a finite dis- 
tance. Shocks occur because the sound speed at the crest of 
the wave is greater than that in the trough. Thus, a wave 
that is launched with a sinusoidal profile eventually steep- 
ens into a saw-tooth form and its energy is dissipated in a 
shock. The number of wavelengths required for the shock 
formation to occur can be ~ c/Ac, where Ac is the differ- 
ence in sound speed between the crest and the trough of 
the wave. This argument ignores various geometrical cor- 
rections that occur in a multi-dimensional disc. In addition, 
it ignores the possibility that dispersive effects could lessen 
the wave steepening (Larson 1990; Papaloizou & Lin 1995). 
This estimate suggests that the propagation distance before 
shocks set in depends on the initial amplitude of the wave, 
i.e. the wave amplitude at the resonance. 



4.2 General wave properties 

Along with the simple two-dimensional wave, a vertically 
isothermal disc admits series of other waves, all of which 
possess vertical motion (Lubow & Pringle 1993). These wave 
modes can be categorised by their vertical velocity structure. 
In a thermally stratified disc, as pointed out by Lin, Pa- 
paloizou & Savonije (1990), the two-dimensional wave does 
not exist because 7p/p is a function of z (cf. equation |^). 
In such a disc, all modes involve vertical motions, as follows 
from Korycansky & Pringle (1995). 

In the spherical geometry of the simulations, z = r cos 9. 
The linear wave modes then have the form 



Q(r,6,t) = Re<^ Q(r,9)exp 



iuit + i / k(r) dr 



(9) 



where Q is any perturbation quantity. The frequency uj is de- 
fined to be positive. The wavenumber k(r) satisfies \kr\ 2> 1 
except near resonances, where this WKB form breaks down. 
The wavenumber can be positive or negative, depending on 
the radial direction (inward or outward) of the phase ve- 
locity. The amplitude Q varies slowly with r. The vertical 
structure of the wave, and the relation between uj and k, 
are determined by an eigenvalue problem, local in r, which 
admits a set of discrete modes. Each mode has a different 
dispersion relation uj = uj(k, r) which must be satisfied at 
every r. 

Using the notation of Ogilvie (1998), these modes can 
be categorised into four classes which are determined by 



the dominant forces involved in their propagation. The f 
(fundamental), p (pressure-dominated), and g (buoyancy- 
dominated) modes propagate where uj > Cl(r), while r modes 
(rotationally dominated modes) propagate where uj < O(r). 
There are only two f modes, f B and f° . Each of the other 
classes contains an infinite sequence of modes with increas- 
ing numbers of nodes in the vertical structure. These may 
be classified as pf , pi, pi, . . . , p?, P2, P§, • ■ • , etc., where 
the number refers to the order of the mode, and 'e' or 'o' 
refers to even or odd symmetry about the mid-plane. 

In this paper, we focus our attention on the following 
modes. 

• the f c mode, which propagates where w > f2(r). In a 
vertically isothermal disc, this mode is the two-dimensional 
mode described in the previous subsection. For other types 
of discs, the f e mode behaves like the two-dimensional wave 
close to the resonance, but behaves like a surface grav- 
ity wave away from the resonance (Ogilvie 1998; Lubow & 
Ogilvie 1998). 

• the Pi mode, which propagates where uj > fi(r)y , 7 + 1. 
Its vertical velocity at resonance is proportional to z (or 
cos 9). 

• the f°/r° mode, which propagates at all radii, but has 
a resonance where uj = CI. Its vertical velocity eigenfunction 
at resonance is independent of z, while its radial velocity is 
proportional to z. The f°/r° mode is a corrugation wave or 
disc warp, where the mid-plane of the disc oscillates up and 
down. 

To excite these modes efficiently, we apply an accelera- 
tion (a r ,ao) (force per unit mass) to our equilibrium discs 
of the following forms. 

(i) a r = A(r) sinujt and ag — for the f° mode, with 
A(r) = Ar" 2 ; 

(ii) a r — and ag — A(r)r cos 9 smujt for the p\ mode, 
with A(r) = Aor~ 3 ; 

(iii) a r — and ag = A(r) sinut for the f°/ r i mode, with 
A(r)=A r~ 2 . 

In each case, the power-law dependence of A on r is chosen 
so that the amplitude of the non-resonant response can be 
reasonably small (linear) at both small and large radii. We 
choose a wave frequency uj = 1 so that the resonance in our 
dimensionless units (uj — CI or uj = + 1) ues at r = 1 

(for the P and corrugation modes) or r = (7 + l) 1 ' 3 « 1.39 
(for the p5 mode). The simulated disc contains these reso- 
nance locations. The acceleration is applied from the begin- 
ning of the calculations and its strength is parameterised by 
Ac 

Each form of the acceleration listed above generally re- 
sults in the excitation of a variety of modes. The total energy 
flux generated at the resonance that is carried by all modes 
for an axisymmetric disc can be determined by adapting 
the torque formula of Goldreich & Tremaine (1979). In Ap- 
pendix B, we derive the total energy flux for each form of 
acceleration listed above. 

In practice, it is more convenient to express the strength 
of the forcing in terms of the spatial peak of the time- 
averaged wave Mach number at the disc mid-plane that oc- 
curs near the resonance, which we denote as M Y . In terms 
of our notation, we have 
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Figure 1. The integrated radial energy flux, |/ r [, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, Ai(r,9 = 
it/2), (bottom) are plotted as functions of radius for the f c mode in a vertically isothermal disc with H/r = 0.1. The results are shown for 
four different wave amplitudes: Air = 0.001 (left), Air = 0.01 (centre-left), Air = 0.03 (centre-right), Ai r = 0.1 (right). Each case was 
performed with different numerical resolutions and/or viscosity to test for convergence: 250 X 111 and g v isc = 4 (long-dashed), 500 X 221 
and <j visc = 4 (triple-dot-dashed), 500 X 221 and <j visc = 0.5 (short-dashed), 1000 X 441 and q visc = 4 (dot-dashed), 1000 X 441 and 
owe = 0.5 (solid). The horizontal dotted line indicates 75 percent of the flux predicted by the Goldreich & Tremaine formula (equation 
pl| ) and the associated Mach number of the wave. The value of 75 percent is that predicted by linear theory to be deposited in the f e 
mode. As can be seen from the figure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude 
cases. In the high-amplitude cases, the radial energy flux is attenuated by dissipation. 



M T = max[M(r,e = w/2)], (10) 

where M is defined by equation (^|) and the maximum 
is taken over a suitable neighbourhood of the resonance 
(0.8 < r < 1.2). This measure of forcing strength is used 
in comparing the simulations involving different disc mod- 
els. In carrying out a simulation with a chosen value of Mi, 
we vary Ao until the desired value of Mr is achieved. 



4.3 Numerical convergence 

As described above, in a vertically isothermal disc, the f c 
mode should propagate outwards from the resonance with- 
out a loss of flux until the wave steepens to form shocks. 
One of the main difficulties in obtaining non-linear results is 
spatially resolving the wave until this wave steepening and 
shock formation occur. If the flow becomes unresolved (i.e. 
it varies on a scale similar to the zone spacing) , the wave will 
dissipate energy. This dissipation occurs in one of two ways. 
In the absence of an artificial viscosity, the energy carried by 
the wave is lost from the calculation. If there is an artificial 
viscosity, as is employed in the calculations presented here, 
the wave energy is converted into thermal energy. Artificial 
viscosity provides a means of simulating shocks (flow discon- 
tinuities) in the code. Numerical dissipation in the ZEUS-2D 
code is proportional to the square of the velocity difference 
across a grid zone. Thus, for an approximately linear wave, 
the numerical dissipation is proportional to q v i BC A 2 /(XN r ), 



where A and A are the amplitude and wavelength of the 
wave, respectively. 

There are several ways to determine whether a wave 
dissipates numerically or physically. First, the resolution of 
the calculations can be increased. If the dissipation is nu- 
merical, the wave will propagate further before it dissipates. 
Secondly, keeping the resolution fixed, the strength of the 
artificial viscosity, g v isc, can be varied. If the dissipation is 
numerical, the wave will propagate further before it dissi- 
pates with lower q v i sc - Note, however, that if g v isc <C 1, 
energy will simply be lost from the calculation instead of 
being converted to thermal energy (e.g. calculations with 
gvisc = 0.1 and g v i sc = 0.01 give almost identical results). 
The third way is to examine the form of the wave to see 
if wave steepening occurs just before the wave dissipates. 
If the wave maintains a sinusoidal profile, the dissipation 
is numerical. In order to determine which of the non-linear 
results are due to physical dissipation, we use all three of 
these indicators. 



5 THE VERTICALLY ISOTHERMAL DISC 

In this section, we study the excitation, propagation and dis- 
sipation of adiabatic (7 = 5/3) waves in a vertically isother- 
mal accretion disc. First, we consider the f c mode (a plane 
wave). Subsequently, we discuss the pf mode and the f°/r° 
(corrugation) mode. 
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Figure 2. The Mach number of the radial motion of the wave 
in the mid-plane of the disc is plotted at a particular instant as a 
function of radius for the f e mode in the vertically isothermal disc 
with H/r = 0.1. Various wave amplitudes are shown: Mr = 0.01 
(top), Mi = 0.03 (middle) and Mr = 0.1 (bottom). For the low- 
amplitude case, Mr = 0.01, the wave maintains its sinusoidal 
profile until numerical dissipation occurs. For the higher ampli- 
tudes, Mr = 0.03 and Mr = 0.1, the wave steepens, shocks and 
dissipates. This indicates that the physical dissipation is numeri- 
cally resolved. 



5.1 f° mode 

5.1.1 Linear theory 

As discussed in Section ^, we adopt a purely radial driving 
force per unit mass of a r (r,t) = A(r) sin uit, in order to effi- 
ciently excite the f c mode. The f c mode corresponds to the 
two-dimensional wave in a vertically isothermal disc. Linear 
theory predicts that the radial forcing of Section ^ should 
lead to the outwardly propagating f c mode receiving 74.5% 
of the total radial energy flux (see Appendix Bl). The in- 
wardly propagating rf mode should receive about 6.0% of 
the total flux. The remainder is goes into g modes (see Ap- 
pendix Bl). 

According to linear theory (Lubow & Pringle 1993), the 
f c mode should occupy the entire vertical thickness of the 
disc and propagate outwards in radius. The mean Mach 
number at the mid-plane of the disc should increase as 
M{r, = 7r/2) oc r 1 / 2 for r S> 1. This weak dependence 
of the Mach number on r means that a low-amplitude wave 
is expected to propagate a large distance before it steep- 
ens into a shock and dissipates. The vertical structure of 
the wave can also be predicted by linear theory. The mean 
Mach number varies vertically as M oc exp(z 2 /5H 2 ), for 
7 = 5/3. Thus, although a wave may be linear near the 



mid-plane, there will be some finite height at which the mo- 
tions are supersonic and we would expect shocks to form. 
However, the majority of the wave energy is contained near 
the mid-plane and therefore this does not diminish the wave 
flux appreciably. 



The r modes will be discussed in detail in Section 5.3 
Briefly, however, the r^ mode should propagate inwards and 
be channelled towards the mid-plane of the disc. 

5.1.2 Low-amplitude, non-linear results 

We turn now to the non-linear hydrodynamic calculations 
performed with ZEUS-2D. Figure |l| plots the magnitude of 
the vertically integrated radial energy flux, J/ r | (equation 
Q), and the mean Mach number (equation q) in the mid- 
plane of the disc, M(r,9 = it/2), as functions of r. They are 
plotted for four values of the Mach number near resonance 
Mr - Notice that the radial energy flux is in fact negative for 
r < 0.7. This effect is due the inwardly propagating r*j; mode 
that is excited along with the f e mode. 

The dotted lines in the upper panels of Figure |l| show 
the contribution from the f c mode to the total radial energy 
flux that is predicted by linear theory. In the lower panels, 
we plot the corresponding Mach numbers that are predicted 
in the mid-plane by linear theory. The non-linear results are 
in excellent agreement with the linear theory. 

Consider the weakest f e mode (left panels, Figure Q). 
Excitation of the f° mode occurs over a finite radial extent 
(a few H) around r — 1. For r > 1, the flux reaches a 
peak value just beyond the resonance. The drop in flux just 
beyond the peak is probably due to the rapid dissipation 
of g modes, which account for about 20% of the total flux. 
The P mode continues to propagate to greater radii and its 
flux remains nearly constant for 2 < r < 8. The value of the 
flux there agrees well with the predicted f° mode flux. The 
dissipation that occurs at large radii is numerical as seen 
by the lack of convergence with increasing resolution. Thus, 
the wave should continue to propagate outwards, beyond the 
grid boundary. 

For r < 1, the r° mode propagates inwards and overlaps 
with the evanescent tail of the P mode. The rf mode is only 
modelled for r ii 0.5 because the edge of the equilibrium disc 
is at Ri — 0.5. The two modes have fluxes of opposite sign, 
and the net flux changes sign near r = 0.7. The peak radial 
energy flux from the simulation in this region is about 3% 
of the total flux. This result is consistent with the expected 
x\ mode negative flux of about 6% of the total overlapping 
with the tail of the positive P mode flux. The r modes will 
be discussed in more detail in section 5. J where an r° mode 
receives 50% of the total flux. 

These results are reinforced by Figure ^, which gives 
contour plots of the wave energy density (equation ||) and 
the energy dissipation rate for the low-amplitude calcula- 
tions (left panels). There is a lot of dissipation at high al- 
titude in the disc in the range 1 < r < 2, as expected for 
the g modes. In the bulk of the disc, there is no significant 
dissipation and the P mode propagates outwards with vir- 
tually no loss of flux (a small amount is lost at high altitudes 
where the motions become supersonic). 

The behaviour of the inwardly propagating r\ mode can 
also be seen in the left-hand panels of Figure M. As predicted 
by linear theory (Lubow & Pringle 1993), the wave is chan- 
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Figure 3. The case considered here is the f e mode in a vertically isothermal disc with H/r = 0.1. The kinetic energy density in the 
wave (logi?, equation ^) is shown in the upper panels. The energy dissipation rate per unit volume is shown in the lower panels as a 
logarithmic grey-scale. The grey-scales cover three orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right 
panels are for amplitude Mr = 0.1. At low amplitude the P mode propagates to large radius throughout the body of the disc unaffected 
by the small amount of dissipation that occurs at high \z\. In the high-amplitude case dissipation occurs in the body of the disc and 
significantly attenuates the wave (see Figure hi). In both cases the grid resolution is N r X Ng = 1000 X 441 and g v isc = 4.0. 



nelled towards the mid-plane. This will be discussed further 



in Section 5.3 



5.1.3 Dependence on wave amplitude 

As the driving amplitude is increased, both numerical and 
physical dissipation occur more quickly. Numerical dissipa- 
tion incr ease s as the square of the amplitude of the wave 
(Section O). Thus, waves with larger values of Mr propa- 



gate a shorter distance before numerical dissipation occurs 
(cf. Figure [j], results with Mr = 0.001, and Mr = 0.01). 

Physical dissipation occurs if the sinusoidal wave steep- 
ens into a shock, which takes place in ~ c/Ac wavelengths 
from the resonance (Section For waves with very small 
values of Mr, steepening requires many wavelengths and 
numerical dissipation sets in first, as we saw in the case for 
Mr = 0.001. Physical dissipation in our simulations was 
obtained for cases with Mr = 0.03 and 0.1. Consider the 
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Figure 4. The integrated radial energy flux, |/ r [, (top) and mean (time-averaged) Mach number at the mid- plane of the disc, M(r, 9 = 
7r/2), (bottom) are plotted as functions of radius for the f c mode in a vertically isothermal disc with H/r = 0.05. The results are 
shown for three different wave amplitudes: Ai T = 0.001 (left), Ai T = 0.01 (centre), Ai T = 0.1 (right). Each case was performed with 
different numerical resolutions to test for convergence: 250 X 111 and <j v i sc = 4 (long-dashed), 500 X 221 and g v isc = 4 (triple-dot-dashed), 
1000 X 441 and q„i y = 4 (dot-dashed). The horizontal dotted line indicates 75 percent of the flux predicted by the Goldreich & Tremainc 
formula (equation Klhand the associated Mach number of the wave. The value of 75 percent is that predicted by linear theory to be 
deposited in the f e mode. As can be seen from the figure, the hydrodynamic results are in excellent agreement with linear theory for the 
lowest amplitude case. In the high-amplitude cases, the radial energy flux is attenuated by dissipation. 



panels for Mr = 0.03 in Figure [j]. Convergence of flux and 
M(r, 9 = 7r/2) is obtained for r 4 with resolutions of 
500 x 221 and 1000 x 441. With an initial amplitude of 
Mr — 0.1, even a resolution of 250 x 111 is sufficient to ob- 
tain physical dissipation. In the Mr = 0.03 case, the wave 
flux has been reduced by an order of magnitude by r ~ 10, 
whereas in the Mr = 0.1 case this level of reduction occurs 
at r w 4. 

To emphasise that this dissipation is physical rather 
than numerical, we plot in Figure ^ the radial velocity 
profile of the wave in the mid-plane for the cases with 
Mr = 0.01,0.03, and 0.1. Steepening of the wave towards 
the profile of a shock-wave is clearly visible in those cases 
with Mr ~ 0.03 and 0.1, whereas in the Mr = 0.01 calcula- 
tion the profile remains sinusoidal indicating only numerical, 
rather than physical, dissipation. This result is consistent 
with Figure g which shows that increasing Mr decreases 
the volume in which the f e mode propagates. 



5.1.4 Dependence on disc thickness 

In Figures ^| and |f| we plot the energy flux and mean Mach 
number M(r, 9 = 7r/2) as functions of r for discs with 
H/r — 0.05 and H/r — 0.2, respectively. The main dif- 
ference between the waves in these discs and those in the 
H/r = 0.1 disc is that the wavelength is proportional to H/r . 
Since the numerical dissipation is tx l/(XN r ) (Section 4.3), 
a wave in a thick disc can be followed to larger radii than in 
a thin disc for the same radial resolution (cf. Mr = 0.01 in 



Figures |lj ^ and q) . Thus , for example , waves with Mr = 0.1 
are well resolved in the hottest disc, H/r = 0.2, even with 
N T x Ng — 250 x 111, and dissipate due to shocks, whereas 
in the coldest disc, H/r = 0.05, even the highest resolutions 
do not give convincing convergence (numerical dissipation 
still plays a significant role). 

Strictly, linear theory assumes H r. However, we find 
that, even with H/r = 0.2, linear theory provides a good de- 
scription of the wave excitation and propagation. A further 
prediction that can be derived from linear theory is that 
the ratio of the Mach number at the resonance to the Mach 
number at r > r ICS should be nearly independent of H/r. 
This is confirmed by the non-linear calculations presented 
here. 



5.2 pi mode 

5.2.1 Linear theory 

As discussed in Section ^j, we adopt a driving force per unit 
mass that is purely in the ^-direction of the form ae(r, 9, t) = 
A{r)r cos 9 sin ut, in order to efficiently excite the p° mode. 
The pf mode is launched at a vertical resonance r ~ 1.39 
and propagates outwards. We predict its radial energy flux 
in Appendix B2. No other wave is expected. 

Unlike the f e mode, the p° mode has vertical structure 
(see Lubow & Pringle 1993). In the vicinity of the resonance, 
the wave energy has a minimum on the mid-plane of the 
disc and two maxima at z/H m 2. This is due to the form of 
the excitation which depends linearly on 2 (Section ^). For 
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Figure 6. The integrated radial energy flux, |/ r |, is plotted as a function of radius for the mode in a vertically isothermal disc with 
H/r = 0.1. The results are shown for three different wave amplitudes: Ai r = 0.001 (left), J\4 r = 0.01 (centre), M T = 0.1 (right). Each 
case was performed with different numerical resolutions to test for convergence: 250 X 111 and g v i sc = 4 Hong-d ashed ), 500 X 221 and 
9visc = 4 (triple-dot-dashed). The horizontal dotted line indicates the flux predicted by linear theory (equation B20). As can be seen 
from the figure, the hydrodynamic results are in excellent agreement with linear theory for the low amplitude cases. At high amplitude, 
the radial energy flux is attenuated by dissipation. 
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Figure 5. The integrated radial energy flux, |/ r |, (top) and 
mean (time-averaged) Mach number at the mid-plane of the disc, 
M(r, 9 = 7r/2), (bottom) are plotted as functions of radius for 
the f e mode in a vertically isothermal disc with H/r = 0.2. The 
results are shown for two different wave amplitudes: Ai r = 0.01 
(left), Air = 0.1 (right). Each case was performed with differ- 
ent numerical resolutions to test for convergence: 250 X 111 and 
Qviec = 4 (long-dashed), 500 X 221 and q v isc = 4 (triple-dot- 
dashed), 1000 X 441 and <7 v j sc = 4 (dot-dashed). The horizon- 
tal dotted line indicates 75 percent of th e flu x predicted by the 
Goldreich & Tremaine formula (equation Bl) and the associated 
Mach number of the wave. The value of 75 percent is that pre- 
dicted by linear theory to be deposited in the f e mode. As can 
be seen from the figure, the hydrodynamic results for the energy 
flux are in good agreement with linear theory for the lowest am- 
plitude case, but for this relatively thick disc the predicted Mach 
number underestimates the computed value by about 30 percent. 
In the high-amplitude case, the radial energy flux is attenuated 
by dissipation. 



r ^ 1.7, this changes to three maxima, one on the mid-plane 
and two at z/H « 3, with two minima at z/H ~ 2. Lubow 
& Pringle (1993) interpreted this as the wave energy rising 
up within the disc as the wave propagates outwards. In fact, 
however, the vertical distribution of the wave energy does 
not change beyond r « 2. The wave propagates outwards 
maintaining the above distribution of energy indefinitely. As 
with the f c mode, the wave is eventually expected to undergo 
steepening into a shock and dissipate, but a low-amplitude 
wave will propagate for a large distance before this takes 
place. 



5.2.2 Low-amplitude, non-linear results 

Figure [] plots the magnitude of the vertically integrated ra- 
dial energy flux, |/ r j (equation as a function of r. We con- 
sider three values of the mean Mach number near resonance 
Mr- As with the case for the f e mode, the low-amplitude pj 
mode propagates outwards from the resonance until most of 
the energy is dissipated numerically. 

A small negative energy flux is seen for r <J 0.9. We 
attribute this to a low-amplitude r mode that is excited at 
the Lindblad resonance at r = 1. In a very thin disc, vertical 
forcing of the type applied here should not excite any modes 
at the Lindblad resonance. However, when H/r — 0.1, there 
is some ambiguity between the meanings of 'horizontal' and 
'vertical' away from the mid-plane. As a result, some launch- 
ing at the Lindblad resonance is possible at the level of 
(H/r) 2 . 

In Figure m we plot the wave energy and energy dissipa- 
tion rate for the low-amplitude pi mode (left panels). The 
vertical distribution of the wave energy density is in excel- 
lent agreement with linear theory in the bulk of the disc with 
only a small amount of dissipation at high altitude where the 
motions become sonic. 



5.2.3 Wave amplitude and disc thickness 

As with the f c mode, the dissipation occurs more rapidly 
with increased driving. Calculations with A4 r = 0.1 resolve 
the wave until physical dissipation due to wave steepening 
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Figure 7. The case considered here is the p| mode in a vertically isothermal disc with H/r = 0.1. The kinetic energy density in the 
wave (log_B, equation is shown in the upper panels covering four orders of magnitude. The energy dissipation rate per unit volume is 
shown in the lower panels as a logarithmic grey-scale covering six orders of magnitude. The left panels are for amplitude _M r = 0.001 
and the right panels are for amplitude M T = 0.1. At low amplitude the p^ mode propagates to large radius throughout the body of 
the disc unaffected by the small amount of dissipation that occurs at high \z\. The p^ mode displays the vertical structure predicted by 
linear theory. In the high-amplitude case dissipation occurs in the body of the disc and significantly attenuates the wave (see Figure |i| . 
In both cases the grid resolution is N r X Ng = 500 X 221 and g v i sc = 4.0. 



occurs, but the lower A4 T cases are still limited by numerical 
dissipation. 

We have not performed calculations of the pi mode 
for discs with different thicknesses since, as we found for 
f e mode, we do not expect any significant dependence of the 
propagation of the pi mode on H/r. As with the f e mode, 
the longer wavelength in hotter discs would make it easier 
to resolve the wave to larger radii. 



5.3 P/r? mode 

5.3.1 Linear theory 

As discussed in Section ^j, we adopt a driving force per unit 
mass that purely in the ^-direction of the form ag(r,t) = 
A(r) sinujt, in order to efficiently excite the f°/ r i mode. 
The f°/r° mode, or corrugation mode, is also launched at 
a vertical resonance although, in a Keplerian disc, this co- 
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Figure 8. The integrated radial energy flux, |/ r |, is plotted as a function of radius for the f°/r° (corrugation) mode in a vertically 
isothermal disc with H/r = 0.1. The results are shown for three different wave amplitudes: JV4 T = 0.001 (left), J\4 T = 0.01 (centre), 
.Mr = 0.1 (right). Each case was performed with different numerical resolutions to test for convergence: 250 X 111 and g v i sc = 4 (long- 
dashed), 500 X 221 and g v i sc = 4 (triple-dot-dashed). The horizontal dotted line indicates the flux predicted by linear theory (equation 
B31). For this case, the inwardly propagating rj mode (r < 1) and the outwardly propagating f° mode (r > 1) carry equal amounts of 
energy. As can be seen from the figure, the hydrodynamic results are in reasonable agreement with linear theory for the low-amplitude 
high-resolution case. At high amplitude, the radial energy flux is attenuated by dissipation. 



incides with the Lindblad resonance making it a hybrid ver- 
tical/Lindblad resonance. In the discs studied here, this hy- 
brid resonance occurs at r = 1. The f°/r° mode is able to 
propagate throughout the entire disc with equal amounts of 
flux in the inwardly propagating rj mode and the outwardly 
propagating f° mode. In Appendix B3, we determine the 
total energy flux produced at this resonance. 

The t° mode is rapidly channelled towards the mid- 
plane of the disc with an associated rapid rise in the Mach 
number of the wave that might be expected to lead to dissi- 
pation of the wave in shocks (Lubow & Pringle 1993). The 
f° mode has two peaks in the vertical distribution of its wave 
energy with a minimum on the disc mid-plane. 



5.3.2 Low-amplitude, non-linear results 

Figure ^ plots the magnitude of the vertically integrated 
radial energy flux, |/ r | (equation Q) as a function of r. We 
consider three values of the mean Mach number near reso- 
nance Air- As predicted by linear theory, this mode propa- 
gates throughout the entire disc. Outward of the resonance 
(r > 1), the wave propagates as an f° mode. Inward of the 
resonance, the wave propagates as an rj mode. An equal 
amount of flux goes into each of the modes, and the mag- 
nitude of the flux is in good agreement with linear theory, 
particularly for the calculations with the highest resolution. 

Once again, for low-amplitude driving, the distance that 
the waves propagate is limited by numerical dissipation 
rather than physical dissipation (shocks). 

The energy of the f° mode is concentrated away from 
the mid-plane of the disc (Figure ^) and its vertical dis- 
tribution is in excellent agreement with linear theory for 
\z/H\ ^ 3. Above this region and particularly within the 
range 1 < r £ 2.5, shocks form and a small amount of 
energy is dissipated. However, the majority of the wave en- 
ergy remains within \z/H\ 3 and propagates outwards 
indefinitely. Lubow & Pringle (1993) noted the rise of wave 
energy away from the mid-plane that occurs from near the 
resonance to r = 3 and interpreted this as the wave energy 
rising into the atmosphere. However, we see here that this 
'rising' does not continue indefinitely; the wave simply oc- 



cupies the bulk of the disc as it propagates to large radii, in 
a similar manner to the f c mode except that it has vertical 
structure. 

As expected from linear theory (Lubow & Pringle 1993), 
the inwardly propagating rj mode is strongly focused to- 
wards the mid-plane of the disc (Figure ^). Both linear the- 
ory and the non-linear results also show that the wavelength 
of the Ti mode decreases towards the centre. 

5.3.3 Wave amplitude 

With strong driving (e.g. A4 T = 0.1), physical dissipation 
through shocks occurs in the simulations and those with 
different resolution converge, at least for the f° mode. By the 
time the f° mode reaches r — 10, the flux has been reduced 
by a factor of ~ 10 3 . The wave is confined to a region very 
near the mid-plane of the disc; shocks are present at large \z\ 
(Figure ^[ right panels) . We have not performed calculations 
of the corrugation mode with other disc thicknesses. 



6 THE POLYTROPIC DISC 

We now turn our attention to polytropic discs. The disc 
studied in this section has an optical depth r = 25000, and 
only a tenuous isothermal atmosphere. Discs that have r = 
5 — 100 (i.e. intermediate between a vertically isothermal 
disc and a very optically thick disc) are briefly described in 
Section 7. 



6.1 P mode 

6.1.1 Linear theory 

As discussed in Section 0, we adopt a purely radial driving 
force per unit mass of a r (r, t) = A(r) sin ait, in order to effi- 
ciently excite the f° mode. The f° mode corresponds to the 
two-dimensional wave in a vertically isothermal disc. Linear 
theory predicts that this radial forcing should lead to the 
outwardly propagating f° mode receiving 98.3% of the to- 
tal radial energy flux (see Appendix Bl). A greater fraction 



12 M. R. Bate et al. 



Figure 9. The case considered here is the f°/r? (corrugation) mode in a vertically isothermal disc with H/r = 0.1. The kinetic energy 
density in the wave (log E, equation ^) is shown in the upper panels covering five orders of magnitude. The energy dissipation rate per 
unit volume is shown in the lower panels as a logarithmic grey-scale covering six orders of magnitude. The left panels are for amplitude 
M T = 0.001 and the right panels are for amplitude Ai T = 0.1. At low amplitude the f° mode propagates outwards from r = 1 to large 
radius throughout the body of the disc unaffected by the small amount of dissipation that occurs at high \z\. The r° mode propagates 
inwards from r = 1 and is channelled towards the mid-plane of the disc. Both the f° mode and the r° mode display the vertical structure 
predicted by linear theory. In the high-amplitude case dissipation occurs in the body of the disc and significantly attenuates the wave 
(see Figure p|). In both cases the grid resolution is N r X Ng = 500 X 221 and q v isc = 4.0. 



goes into the f c mode in the polytropic disc than in the ver- 
tically isothermal disc. The remainder goes into an inwardly 
propagating rl mode. 

The linear behaviour of the f e mode in a polytropic disc 
has been studied by Korycansky & Pringle (1995), Lubow & 
Ogilvie (1998) and Ogilvie & Lubow (1999). Although the f e 
mode is excited throughout the entire vertical extent of the 



disc at the resonance, it is channelled towards the surface of 
the polytropic disc, unlike in the vertically isothermal disc. 
The rate of channelling (the distance from the resonance at 
which the wave energy becomes concentrated away from the 
disc mid-plane) depends only on the degree of vertical strat- 
ification, which is determined by the index of the polytrope, 
and the azimuthal wavenumber m, which is fixed at m = 
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for the axisymmetric case. The rate is independent of the 
disc thickness H, in the limit that H <C r. This concentra- 
tion of the wave energy into a small volume of cool, low- 
density gas near the disc surface (or atmosphere) leads to a 
rapid increase of the Mach number of the wave with radius. 
Thus we expect that waves in polytropic discs will undergo 
shocks more quickly (closer to the resonance) than waves 
in the isothermal case described earlier. Lubow & Ogilvie 
(1998) speculated that once the Mach number of the wave 
became of order unity, the wave would shock and thus, such 
waves would not propagate over large distances. One of the 
aims of the non-linear hydrodynamic calculations is to de- 
termine where this dissipation occurs. 

6.1.2 Low-amplitude, non-linear results 

We turn now to the non-linear hydrodynamic calculations 
performed with ZEUS-2D. Figure hd plots the magnitude of 
the vertically integrated radial energy flux, |/ r | (equation 
as a function of r. We consider four values of the Mach num- 
ber near resonance M r . The calculations with Mi — 0.001 
do not have sufficient resolution to follow the f e mode (r > 1) 
until physical dissipation occurs. This can be seen by the 
absence of convergence of the flux in Figure [HJ Because 
of the higher concentration of wave energy compared with 
the vertically isothermal case, it is more difficult to pro- 
vide adequate resolution. However, it is clear that the wave 
propagates to at least r = 2.3 with little dissipation, over a 
distance that is many times greater than the disc thickness 
H. Notice that the radial energy flux is negative for r < 0.7. 
This effect is due the inwardly propagating rf mode that is 
excited along with the f e mode. As predicted by linear the- 
ory, most of the flux goes into an outwardly propagating f c 
mode. The magnitude of this flux is indistinguishable from 
the linear prediction. We find w 1.2% goes into an inwardly 
propagating rf mode (Figure [l^, Mi = 0.001), similar to 
the linear prediction. 

The vertical thermal stratification of the polytropic disc 
has a strong effect on the distribution of wave energy of the 
f c and rf modes compared to the isothermal disc (Figure 
[ll], top left). As predicted by linear theory, instead of occu- 
pying the entire thickness of the disc, as in the isothermal 
case, the energy in the f c mode is channelled towards the 
surface of the disc. It effectively becomes a wave that trav- 
els along the surface of the disc. For the rf mode, rather 
than being channelled towards the mid-plane, it continues 
to occupy the entire thickness of the disc as it propagates 
towards the central object. This behaviour differs from the 
vertically isothermal case, where the rf mode is channelled 
towards the mid-plane. The reason for this difference is that 
the channelling of the rf mode is due to vertical buoyancy 
(Lubow & Pringle 1993). The difference in behaviours is due 
to the vanishing buoyancy frequency throughout the verti- 
cally polytropic disc, unlike the vertically isothermal case. 

The wave energy and dissipation for our weakest driv- 
ing calculation, Mi — 0.001, are shown in Figure [ll] (left 
panels). The dissipation is only suggestive, since numerical 
convergence has not occurred at larger radii. However, the 
channelling of the f e mode to the surface of the disc is obvi- 
ous. The wave dissipates at r « 2.3 in the highest resolution 
calculation (N r x Ng = 1156 x 365) and at r w 2.2 with half 
of this resolution (N r x N = 1000 x 183). 
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Table 1. The radius r at which the f° mode dissipates is tab- 
ulated for polytropic discs in which the calculations resolve the 
shock dissipation. Values are given for different disc thicknesses, 
H/r, and for different initial wave amplitudes M r . The '+' in- 
dicates that this value is a lower limit. Note that the radius of 
dissipation is almost independent of the disc thickness as pre- 
dicted by the theory of wave channelling. 

The low-amplitude rf mode propagates towards the cen- 
tre of the disc (Figure [111 (left panels). The zig-zag pattern 
that can be seen in the distribution of wave energy is due 
to the superposition of the rf mode (for which the wave 
energy has a minimum on the mid-plane) and the evanes- 
cent tail of the f e mode. The wavelength of the rf mode is 
equal to the wavelength of this zig-zag pattern and decreases 
rapidly with decreasing radius. Linear theory predicts that 
the wavelength should be oc r 5//2 and this is confirmed by the 
non-linear calculations. As the wave propagates towards the 
centre of the disc it dissipates numerically because the wave 
becomes unresolved (the size of the grid zones decreases log- 
arithmically, while the wavelength decreases more rapidly). 

6.1.3 Dependence on wave amplitude 

As described above, the resolution of the calculations with 
Mi = 0.001 is not sufficient to resolve physical dissipation 
of the f° mode due to shocks. However, calculations with 
Mi > 0.01 are able to resolve the f° mode until physical 
dissipation occurs. Considering the convergence of |/ r | in 
Figure |Io[ we find that physical dissipation is marginally 
resolved in the Mi = 0.01 calculations with N r x Ng — 
1156 x 365 (there is very little difference in the radius at 
which |/ r | rapidly decreases between this calculation and 
that with N r x Ng = 1000 x 183). Calculations with Mi > 
0.03 are easily resolved with N r xN g = 1000 x 183. The radii 
to which these waves propagate before physical dissipation 
occurs are given in Table hi 

Compared with the case of the vertically isothermal 
disc, the wave propagates a shorter distance before dissi- 
pating, although still much greater than H . The increased 
non-linearity is due to the increase in the magnitude of ve- 
locity perturbations and the decrease in local sound speed 
that occurs as the wave is channelled into the small region 
near the disc surface. 

In Figure we plot the height above the mid-plane of 
the peak of the wave energy density as a function of radius 
and compare it with the prediction from linear theory. We 
also plot the mean Mach number of the wave along this peak 
in wave energy M (r, 6) and compare this to linear theory. 
Away from the resonance (r <; 1.3), where the linear theory 
should give an accurate description, the agreement between 
the linear theory and the non-linear hydrodynamic calcu- 
lations is excellent for both the location of the peak in the 
wave energy density and the Mach number of the wave along 
this maximum. It appears that shock dissipation occurs once 
the mean Mach number M exceeds approximately 0.15 (i.e. 
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Figure 10. The integrated radial energy flux, |/ r |, is plotted as a function of radius for the f e mode in a polytropic disc with H/r = 0.1 
and t = 25000. The results are shown for four different wave amplitudes: Air = 0.001 (left), Air = 0.01 (centre-left), Ai T = 0.03 
(centre-right), Air = 0.1 (right). Each case was performed with different numerical resolutions to test for convergence: 250 X 45 and 
9visc = 4 (long-dashed), 500 X 91 and g v i sc = 4 (triple-dot-dashed), 1000 X 183 and g v i sc = 4 (dot-dashed), 1156 X 365 and a..;,— = 4 
(solid). The horizontal dotted line indicates 98 percent of the flux predicted by the Goldreich & Tremaine formula (equation Bll). The 
value of 98 percent is that predicted by linear theory to be deposited in the f e mode. As can be seen from the figure, the hydrodynamic 
results are in excellent agreement with linear theory for the low-amplitude cases. In the high-amplitude cases, the radial energy flux is 
attenuated by dissipation. 



Figure 11. The case considered here is the f e mode in a polytropic disc with H/r = 0.1 and r = 25000. The kinetic energy density 
in the wave (logi?, equation ^) is shown in the upper panels. The energy dissipation rate per unit volume is shown in the lower panels 
as a logarithmic grey-scale. The grey-scales cover four orders of magnitude. The left panels are for amplitude Air = 0.001 and the right 
panels are for amplitude Air = 0.1. At low amplitude the f° mode propagates outwards from r = 1 and becomes increasingly channelled 
to z H. The dissipation occurs at r a 2.2 as the wave becomes non-linear (upper panel, Figure |l3[). In the high-amplitude case, wave 
channelling of the f e mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern 
visible at r < 1 is generated by a superposition of the evanescent tail of the f e mode and a low-amplitude inwardly propagating r mode. 
In both cases the grid resolution is AV X Ng = 1000 X 183 and <? v isc = 4.0. 
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Figure 12. Comparison of the non-linear P mode with linear 
theory for the polytropic disc with H/r = 0.1 and r = 25000. 
Top panel: We plot as a function of radius, the height, z, in the 
disc where the wave energy density reaches a maximum. The lin- 
ear prediction is given by the dotted line and the low-amplitude 
Mr = 0.001 computation is given by the solid line. For higher 
amplitude waves the maxima lie roughly along the same curve, 
but attempting to extract them from the calculations results in 
curves that are so noisy that plotting them makes the figure un- 
intelligible. Bottom panel: The time-averaged wave Mach number 
at the location of the peak in wave energy density (as defined in 
the upper panel) is plotted as a function of radius. Four cases are 
shown: non-linear calculations with (lower to upper solid lines) 
Mr = 0.001, Mt = 0.01, Mr = 0.03, Mr = 0.1 and the cor- 
responding linear predictions. The numerical results display ex- 
cellent agreement with linear theory both for the location of the 
maximum wave energy and for the Mach number of the wave at 
these maxima, until the wave dissipates in shocks (Mr = 0.01, 
Mr = 0.03, Mi = 0.1) or numerical dissipation (Mt = 0.001). 



the peak-to-trough mean Mach number exceeds « 0.4). Note 
that the Mach number reached by the Mr = 0.001 calcula- 
tion at r — 2.5 is close to this value, implying that even the 
lowest-amplitude wave may be marginally resolved with the 
resolution of N r X Ng = 1156 X 365. 

Finally, we note that the efficiency of the dissipation de- 
pends on Mr, in that stronger waves result in significant flux 
beyond the dissipation radius (cf. Mr = 0.01 and Mr = 0.1 
m Figure [u]). This might be because the strong wave dissi- 
pation modifies the equilibrium disc structure in a manner 
that makes non-linear dissipation more difficult. 
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Figure 13. The Mach number of the radial motion of the wave at 
the height in the disc at which the wave energy density peaks (see 
Figure |l2[ ) is plotted at a particular instant as a function of radius 
for the f 3 mode in the polytropic disc with H/r = 0.1 and r = 
25000. Various wave amplitudes are shown: Mr = 0.001 (top), 
Mr = 0.01 (middle) and Mr = 0.03 (bottom). For the lowest 
amplitude case, Mr = 0.001, the wave maintains its sinusoidal 
profile, but with increasing amplitude as it is channelled towards 
the surface, until numerical dissipation occurs. For the higher 
amplitudes, Mr = 0.01 and Mr = 0.03 the wave becomes non- 
linear very rapidly, shocks and dissipates. 



6.1.4 Dependence on disc thickness 

We have seen that in a polytropic disc the f e mode is 
channelled towards the surface of the disc and dissipates 
due to wave steepening when the mean Mach number at 
the location where the wave energy is concentrated exceeds 
M ~ 0.15. We examine here how this behaviour depends 
on the thickness of the disc. Recall also that that the linear 
theory has been developed under a thin-disc approximation. 

Figures [l4] and [IB] correspond to Figure but for val- 
ues of H/r equal to 0.05 and 0.2 respectively. In Table |l], we 
give the radius of dissipation of the f c mode for those cal- 
culations in which the wave dissipates due to shocks (rather 
than numerically). Although H/r is varied by a factor of 4, 
the radius at which the wave dissipates is almost indepen- 
dent of H/r. This result is consistent with the expectations 
of wave channelling, since channelling occurs over a distance 
that is independent of H/r (Lubow & Ogilvie 1998). 

The main difference between discs of different thickness 
is that thinner discs seem to result in more efficient dissipa- 
tion of the wave. For example, consider the calculations with 
Mr = 0.01. The radial energy flux drops by nearly a factor 
of ~ 1000 between r = 1.7 and r = 2.5 for the H/r = 0.05 
disc, whereas it only drops by a factor of » 10 for H/r — 0.1 
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Figure 14. The integrated radial energy flux, |/ r |, is plotted as a function of radius for the f c mode in a polytropic disc with H/r = 0.05 
and t = 25000. The results are shown for four different wave amplitudes: Air = 0.001 (left), Air = 0.01 (centre-left), A4 T = 0.03 
(centre-right), Air = 0.1 (right). Each case was performed with different numerical resolutions to test for convergence: 250 X 45 and 
9visc = 4 (long-dashed), 500 X 91 and g v i sc = 4 (triple-dot-dashed), 1000 X 183 and g v i sc = 4 (dot-dashed), 1156 X 365 and ffluipc = 4 
(solid). The horizontal dotted line indicates 98 percent of the flux predicted by the Goldreich & Tremaine formula (equation Bll). The 
value of 98 percent is that predicted by linear theory to be deposited in the f e mode. As can be seen from the figure, the hydrodynamic 
results are in excellent agreement with linear theory for the low-amplitude cases. In the high-amplitude cases, the radial energy flux is 
attenuated by dissipation. 



Amplitude, m r =0.01 



Amplitude, T!l=0.03 



Amplitude, 711=0.] 






10 " 


3 




L_ 


10" 6 


CD 




3) 

l5 




O 


10" 7 


*~0 
D 
Ct 




~o 


10" 8 


O 




:ntegi 


10-9 



0.0 0.5 1.0 1.5 2.0 Z.5 3.1 
Radius, r 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 
Radius, r 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
Radius, r 



Figure 15. The integrated radial energy flux, |/ r |, is plotted as a function of radius for the f e mode in a polytropic disc with H/r = 0.2 
and t = 25000. The results are shown for three different wave amplitudes: .Mr = 0.01 (left), .Mr = 0.03 (centre), .Mr = 0.1 (right). 
Each case was performed with different numerical resolutions to test for convergence: 250 X 45 and g v isc = 4 (long-dashed), 500 X 91 
and g v isc = 4 (triple-dot-dashed), 1000 X 183 and g v isc = 4 (dot-dashed), 1156 X 365 an d q, ,;„. = 4 (solid). The horizontal dotted line 
indicates 98 percent of the flux predicted by the Goldreich & Tremaine formula (equation |Bjj| ). The value of 98 percent is that predicted 
by linear theory to be deposited in the f e mode. As can be seen from the figure, the hydrodynamic results are in reasonable agreement 
with linear theory for the most highly resolved low-amplitude case. In the high-amplitude cases, the radial energy flux is attenuated by 
dissipation. 



and only a factor of ~ 3 for H/r = 0.2. The reason for this 
dependence of the dissipation efficiency on the disc thickness 
is not clear. 

In summary, linear theory, although derived for thin 
discs, gives a good description of the wave propagation even 
in discs as hot as H/r = 0.2, although the dissipation of the 
waves may be less efficient for thicker discs. 



6.2 p c L mode 

6.2.1 Linear theory 

As discussed in Section ^j, we adopt a driving force per unit 
mass that is purely in the ^-direction of the form ag (r, 9,t) = 
A(r)r cos sin uit, in order to efficiently excite the Pi mode. 
As described earlier, the Pi mode is launched at a vertical 
resonance at r = 1.39 and propagates outwards. No other 
wave is expected and the radial energy flux of the pi mode is 
predicted in Appendix B2. In a vertically polytropic disc, the 
Pi mode energy is expected to channel to the disc surface. 



6.2.2 Low-amplitude, non-linear results 

Figure |l^ plots the magnitude of the vertically integrated 
radial energy flux, |/ r | (equation Q) as a function of r. We 
consider three values of the mean Mach number near res- 
onance Air- The wave excitation is similar to that in the 
isothermal disc with most of the energy going into an out- 
wardly propagating p\ mode which is excited at r = 1.39 
and a small fraction producing an inwardly propagating r x 
mode which is excited at r = 1. We find that the ratio of 
the two fluxes is w 1.1% (Figure [k| Mr = 0.01). As noted 
previously, this effect is of order (H/r) 2 and should vanish 
in the limit of a thin disc. 

In the isothermal disc, the peak of wave energy of the 
Pi mode in the vicinity of the resonance (r = 1.39) occurs 
off the mid-plane for 1.0 & r 1.7. (There is one maxi- 
mum above, and one below the mid-plane). This structure 
changes as the wave propagates to larger radii to a peak on 
the mid-plane and two off-mid-plane peaks for r <; 1.7. In 
the polytropic disc (Figure [n], top left), this general struc- 
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Figure 16. The integrated radial energy flux, |/ r [, is plotted as a function of radius for the p^j mode in a polytropic disc with H/r = 0.1 
and r = 25000. The results are shown for three different wave amplitudes: M T = 0.001 (left), Ai r = 0.01 (centre), M r = 0.1 (right). 
Each case was performed with different numerical resolutions to test for convergence: 250 X 45 and <j v isc = 4 (long-dashed), 500 X 91 and 
9visc = 4 (triple- dot-d ashed) . 1000 X 183 and q v i sc = 4 (dot-dashed). The horizontal dotted line indicates the flux predicted by linear 
theory (equation B20). As can be seen from the figure, the hydrodynamic results are in excellent agreement with linear theory for the 
low-amplitude cases. At high amplitude, the radial energy flux is attenuated by dissipation. 



Figure 17. The case considered here is the p^j mode in a polytropic disc with H/r = 0.1 and r = 25000. The kinetic energy density 
in the wave (log E, equation ^) is shown in the upper panels covering three orders of magnitude. The energy dissipation rate per unit 
volume is shown in the lower panels as a logarithmic grey-scale covering four orders of magnitude. The left panels are for amplitude 
M T = 0.001 and the right panels are for amplitude M r = 0.1. At low amplitude the p^ mode propagates outwards from r ~ 1.39 and 
becomes increasingly channelled to z H. The mode displays vertical structure as predicted by linear theory (cf. the vertically isothermal 
disc, Figure The dissipation occurs at r a 2.5 as the wave becomes non-linear. In the high-amplitude case, wave channelling of the 
p® mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern visible at r < 1 is 
generated by a superposition of the evanescent tail of the f e mode and a low-amplitude inwardly propagating r mode. In both cases the 
grid resolution is JV,- X Ng = 1000 X 183 and q v i sc = 4.0. 
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Figure 18. The integrated radial energy flux, |/ r |, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, 
_M(r, 9 = tt/2), (bottom) are plotted as functions of radius for the f°/r° (corrugation) mode in a polytropic disc with H/r = 0.1 and 
t = 25000. The results are shown for three different wave amplitudes: M r = 0.001 (left), M T = 0.01 (centre), M r = 0.1 (right). Each 
case was performed with different numerical resolutions to test for convergence: 250 X 45 and q vlac = 4 (long-dashed), 500 X 91 and 
9visc = 4 (triple- dot-d ashed) , 1000 X 183 and q v i sc = 4 (dot-dashed). The horizontal dotted line indicates the flux predicted by linear 
theory (equation B3f). For this case, the inwardly propagating r° mode (r < 1) and the outwardly propagating f° mode (r > 1) carry 
equal amounts of energy. As can be seen from the figure, the hydrodynamic results are in reasonable agreement with linear theory for 
the low-amplitude high-resolution case. At high amplitude, the radial energy flux is attenuated by dissipation. It can be seen that the 
Mach number of the rj mode increases as it propagates inwards caused by the rapid decrease of the group velocity (or radial wavelength, 
Figure Fi). 



ture is repeated, but with the added complexity of the wave 
simultaneously being channelled to the disc surface in a sim- 
ilar manner to the f e mode. Thus, in the vicinity of the res- 
onance (1.0 ^ r 1.7) the wave energy has a minimum 
on the mid-plane as in the isothermal disc. At r <; 1.7, the 
two off-mid-plane peaks which were present in the isother- 
mal disc are near the surface of the polytropic disc. The 
peak which was on the mid-plane in the isothermal disc is 
also present in the polytropic disc for 1.7 ^ r ^ 2.1, but is 
rapidly channelled towards the surface of the disc for r ^ 2.1 
in a manner which is similar to the channelling of the f e 
mode. For r <; 2.5, there is very little wave energy near the 
mid-plane of the disc. 

The rf mode, which is excited at r = 1 and propagates 
inwards, behaves in the same way as the rf mode which 
was excited by the forcing for the f c mode. It occupies the 
entire thickness of the disc and shows no sign of physical 
dissipation as it propagates to the centre. It is dissipated 
numerically because of the rapid decrease of wavelength. 



6.2.3 Wave amplitude 

We only have sufficient resolution to find the radius at which 
physical dissipation occurs for the case with Mr = 0.1 (Fig- 
ure |l^). However, in cases with the same radial resolutions, 
the radius of dissipation is larger for the pi mode than for 



the f c mode (c.f. Figure ^) because the wave is excited at 
r — 1.39 instead of r — 1. 



6.3 P/r? mode 

6.3.1 Linear theory 

As discussed in Section ^j, we adopt a driving force per unit 
mass that is purely in the ^-direction of the form ag(r,t) = 
A(r) sin wt, in order to efficiently excite the f°/r° mode. The 



general properties of the waves were described in Section 5.3 
for the vertically isothermal case. In the polytropic case, the 
primary new feature is that wave channelling occurs. 



6.3.2 Low-amplitude, non-linear results 

Figure ^ plots the magnitude of the vertically integrated 
radial energy flux, |/ r | (equation Q), and the mean Mach 
number (equation ^) in the mid-plane of the disc, M(r, 6 = 
7r/2), as functions of r. They are plotted for three values 
of the mean Mach number near resonance A4 V . The wave 
excitation is similar to that in the isothermal disc with 50% 
of the flux going into each of the inwardly and outwardly 
propagating modes. 

In the isothermal disc, the wave energy of the f° mode 
peaked off the mid-plane. In the polytropic disc (Figure |l9[ 
top left), this is apparent near the resonance (1 ^ r ^ 1.3), 
but for r <; 1.3, channelling of the wave to the disc surface 
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Figure 19. The case considered here is the f°/r° (corrugation) mode in a polytropic disc with H/r = 0.1 and r = 25000. The kinetic 
energy density in the wave (log_E, equation ^) is shown in the upper panels. The energy dissipation rate per unit volume is shown in the 
lower panels as a logarithmic grey-scale. The grey-scales cover four orders of magnitude. The left panels are for amplitude _M r = 0.001 and 
the right panels are for amplitude _M r = 0.1. At low amplitude the f° mode propagates outwards from r = 1 and becomes increasingly 
channelled to z ~ H. Dissipation of the f° mode occurs at r r; 2.2 as the wave becomes non-linear. In the high-amplitude case, wave 
channelling of the f° mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern 
visible at r < 1 is generated by a superposition of the evanescent tail of the f° mode and the r° mode that propagates inwards from 
r = 1. Note that as the r° mode propagates inwards its radial wavelength (and hence its group velocity) decreases. In both cases the 
grid resolution is N r X No = 1000 X 183 and q visc = 4.0. 



dominates the distribution of wave energy. In fact, apart 
from the minimum in the wave energy on the mid-plane 
for 1 ^ r ^ 1.3, the distribution of wave energy is almost 
identical to that of the f c mode (cf. Figure [ll], top left). 

The r° mode, which is excited at r = 1 and propagates 
inwards, occupies the entire thickness of the disc as did the 
rf mode which was excited by the forcing for the f c and 
pj modes. The mean Mach number at mid-plane increases 
rapidly with decreasing radius for the rj mode (Figure |l8[ 
lower panels). In the lowest amplitude, highest resolution 
calculation, it increases by a factor of ~ 4 in going from the 
resonant radius r = 1 to radius r = 0.3. In reality, this would 
result in the formation of shocks and dissipation of the wave 
as it approaches the centre of the disc. However, in the cal- 
culations presented here, we do not have sufficient resolution 
to resolve physical dissipation as can be seen by the lack of 
convergence of |/ r | with increasing resolution (Figure |is| ) . 
Instead, the r° mode is dissipated numerically because of 
the rapid decrease of wavelength. 



6.3.3 Wave amplitude 

For a given wave amplitude, it appears that the f° mode 
dissipates at about the same radius as the f c mode. This 
confirms what was found in comparing the distance of prop- 
agation of the f G and pf modes, namely that the channelling 
of the waves towards the disc surface and their resulting ra- 
dius of dissipation occurs at about the same distance from 
resonance. This result is probably a consequence of the fact 
that wave channelling occurs when kH is greater than unity 
(Lubow & Ogilvie 1998). This condition is determined by 
the wave dispersion relations and occurs at roughly the same 
distance from resonance (namely ~ r) for the modes consid- 
ered here. 



7 DISCS WITH t = 5 - 100 

Thus far, we have considered vertically isothermal discs, 
with optical depth r ~ 0, or polytropic discs with tenuous 
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Figure 20. We consider the f e mode in polytropic discs with H/r = 0.1 and varying optical depth. From top to bottom: r = 5, r = 10, 
r = 30, and r = 100. In each case, the amplitude is Ai T = 0.01. The kinetic energy density in the wave (log E, equation 5) is shown in the 
left-hand panels. The energy dissipation rate per unit volume is shown in the middle panels as a logarithmic grey-scale. The grey-scales 
cover four orders of magnitude. The right panels give the integrated radial energy flux as functions of radius. The horizontal dotted lines 
indicate the flux predicted by linear theory for the f e mode. As can be seen from the figure, the predicted fluxes for discs of intermediate 
optical depth are in excellent agreement with linear theory. As the optical depth increases, there is a gradual change from the vertically 
isothermal behaviour (i.e. the wave occupying the entire thickness of the disc, top) to wave channelling (bottom). This is predicted by 
linear theory (Ogilvie &; Lubow 1999). As wave channelling becomes more important, the attenuation of the wave increases. In each case 
the grid resolution is N r X Ng = 500 X 183 and the disc is modelled over 6 vertical scale- heights. 



isothermal atmospheres and r ~ 25000. Wave propagation 
and dissipation are quite different for these two types of 
discs. In this section we briefly examine cases of intermedi- 
ate optical depth, as denned by equation |^. This has been 
studied in linear theory by Ogilvie & Lubow (1999). They 
found that as the optical depth of the disc is decreased, out- 
wardly propagating waves are channelled into thicker layers. 
This results in a milder increase of the Mach number as the 
wave propagates outwards. 

In this section, we consider only the f e mode. We per- 
form calculations with discs that have r = 5, 10, 30 and 100. 



Their individual parameters were chosen so that the thick- 
ness of the polytropic layer relative to r is the same for each 
disc {H/r £s 0.1). Each calculation had Mi = 0.01. 

The results are provided in Figure |2^. As predicted by 
linear theory, there is a gradual change from the vertically 
isothermal behaviour (i.e. the wave occupying the entire 
thickness of the disc, top) to wave channelling (bottom) as 
the optical depth is increased. For discs of increasing optical 
depth, a greater fraction of the f c mode's energy is chan- 
nelled to the surface of the disc where it dissipates. Thus, 
the fraction of the f° mode flux that propagates past r ~ 2 
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decreases (right panels). The flux in the r = 10 disc is only 
reduced by a factor of « 2, while the f c mode loses more 
than 80% of its flux if r <; 100. These results demonstrate 
that wave channelling is important even in discs of quite low 
optical depth. 

8 CONCLUSIONS AND DISCUSSION 

8.1 Comparison with linear theory 

We have performed non-linear numerical simulations of the 
excitation, propagation and dissipation of axisymmetric hy- 
drodynamic waves in an accretion disc. The disc is Keple- 
rian, vertically resolved and inviscid apart from an artificial 
numerical bulk viscosity that is required in order to resolve 
shocks. The waves are excited by applying a temporally pe- 
riodic acceleration throughout the computational domain, 
and are launched resonantly at radii where the forcing fre- 
quency coincides with a natural oscillation frequency of the 
disc. We have examined Lindblad resonances, vertical res- 
onances and hybrid vertical/Lindblad resonances. In each 
case the waves propagate radially through the disc away 
from the resonance, changing continuously in vertical struc- 
ture as they do so, and ultimately dissipate. 

This mechanism of generating the waves is directly com- 
parable with the way that non-axisymmetric waves are ex- 
cited in discs by tidal forcing from an orbiting companion in 
a binary or protoplanetary system. A low-mass companion 
such as a sub- Jovian planet orbiting a star exerts its tidal 
influence on the disc mainly through launching f c -mode 
waves at the Lindblad resonances (Goldreich & Tremaine 
1980). Lindblad resonances are also important in systems of 
larger mass ratio, such as circular orbit binary star systems 
of extreme mass ratio (Lin & Papaloizou 1979) and eccentric 
orbit binaries (Artymowicz & Lubow 1994). Related reso- 
nances play a role in the growth of the eccentric disc mode 
in superhump binaries (Lubow 1991) and possibly binaries 
with brown dwarf secondaries (Papaloizou, Nelson & Mas- 
set 2001). Vertical resonances are an important aspect of the 
tidal interaction in close binary stars (Lubow 1981; Ogilvie 
2002), where the Lindblad resonances are all excluded from 
the disc. The corrugation-mode resonances we have studied 
are the axisymmetric equivalent of the bending-mode res- 
onances observed in Saturn's rings (Shu, Cuzzi & Lissauer 
1983). 

Wherever the waves are linear, our results are in ex- 
cellent and detailed agreement with the linear theory de- 
veloped by Lubow & Pringle (1993), Korycansky & Pringle 
(1995), Ogilvie (1998), Lubow & Ogilvie (1998) and Ogilvie 
& Lubow (1999). Elements of this theory can also be found 
in the work of Loska (1986) and Okazaki, Kato & Fukue 
(1987). The disc acts as a waveguide that allows a variety of 
wave modes to propagate in the radial direction. The modes 
can be classified similarly to stellar oscillations, and we have 
discussed all four classes (f, p, r and g) to some extent in this 
paper. Each mode has a dispersion relation that connects the 
wave frequency and the radial wavenumber at every radius. 
In general, waves of a given frequency can propagate only in 
restricted radial intervals, bounded by turning points where 
the wavenumber vanishes. The turning points correspond to 
the resonant radii where waves are launched by an applied 
periodic forcing. 



The linear theory is confirmed by the simulations in a 
number of respects. First, the concept of the disc as a wave- 
guide that supports discrete propagating modes is validated. 
The waves are vertically confined by the continuous strat- 
ification of the disc, not by zero-density surfaces, artificial 
boundaries or other discontinuities. The solutions exhibit 
the same 'separation of variables' that also features in the 
linear analysis. That is, the waveform is an oscillatory func- 
tion of radius (and time) multiplied by a vertical profile 
that changes slowly with radius. The fact that the veloc- 
ity amplitudes of the linear solutions may formally diverge 
high in the atmosphere does not invalidate the rest of the 
solution when the energy density of the wave is vertically 
confined in the bulk of the disc. Secondly, the energy fluxes 
imparted to the different wave modes through launching at 
different types of resonance are in excellent agreement with 
the predictions of linear theory presented in Appendix B. 
This analysis of resonant wave launching is closely related 
to the well-known resonant torque formula of Goldreich & 
Tremaine (1979) and to the analysis of vertical resonances 
by Lubow (1981). Finally, the detailed waveforms predicted 
by linear theory are confirmed in the simulations (see, e.g., 
Figure ^l] and the next subsection), even when the disc is 
not especially thin. 

8.2 Wave channelling versus refraction 

The continuous change of the vertical profile of a wave as 
it propagates radially is a process we have termed 'wave 
channelling' (Lubow & Ogilvie 1998). The energy density 
of the wave can be channelled either towards the surfaces of 
the disc or towards the mid-plane as it propagates away from 
its site of launching. The results of this paper provide ample 
evidence for both types of wave channelling. In particular, 
the f and p modes in a disc with a vertical temperature 
gradient are channelled towards the surfaces (e.g. Figures 
0,0 and [is]), while the r modes in a disc with a vertical 
entropy gradient are channelled towards the mid-plane (e.g. 
Figures g and|^) as predicted by Lubow & Pringle (1993). 

The wave energy of the f e mode, which is the principal 
mode launched at a Lindblad resonance, rises towards the 
surfaces of a thermally stratified disc. Close to the resonance, 
the P mode occupies the full vertical extent of the disc and 
behaves compressibly. As the mode propagates away from 
the resonance, its energy within the disc midplane region 
rises and becomes confined (or channelled) to a layer at the 
base of the disc atmosphere. The extent of channelling de- 
pends on the degree of thermal stratification. In this regime, 
the wave becomes a surface gravity wave and behaves incom- 
pressibly (Ogilvie 1998; Lubow & Ogilvie 1998). The process 
by which this occurs is wave channelling and not acoustic re- 
fraction. This is demonstrated as follows. 

If refraction were responsible (e.g. Lin et al. 1990a,b), 
the wave would be directed upwards into the isothermal at- 
mosphere after travelling a radial distance of order H. This 
is not observed in the simulations. Plots of the wave energy 
for low-amplitude waves in polytropic discs (Figures jn], ^ 
and p^| ) clearly show that the waves propagate along the 
base of the atmosphere before dissipating in shocks. Propa- 
gation does not switch from horizontal to vertical, as would 
be expected if simple refraction were involved. In fact, the 
f e mode is launched even in an incompressible fluid with a 
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Figure 21. Comparison of the simulations with linear theory, and comparison of wave channelling with acoustic refraction. In the upper 
grey-scale figure, we plot a snapshot of the radial velocity multiplied by the square root of the density, v T% /~p, from a simulation involving 
an f e mode in a polytropic disc that joins smoothly on to a low-mass isothermal atmosphere (H/r = 0.1, r = 25000). (The square of 
this quantity is the radial part of the instantaneous wave energy density, c.f. equation 5). The amplitude of the wave is A4 r = 0.001 
and the grid resolution is N r X Ng = 1000 X 183. In the lower grey-scale figure, we plot the corresponding semi-analytical prediction of 
the linear theory of wave channelling, ignoring the r mode and using the WKB approximation in the radial direction (which results in 
a mild singularity at the Lindblad resonance). The agreement is excellent until dissipation occurs. Over the upper grey-scale, we trace 
the rays that would result from the refraction of an initially vertical wavefront at the resonant radius. The refracted wavefronts, which 
are orthogonal to the rays, would rapidly become severely tilted and the wave would be predicted to propagate almost vertically into 
the isothermal atmosphere. However, the simulation shows that the wavefronts in fact remain nearly vertical as the wave propagates 
horizontally along the base of the disc atmosphere (shown by the dashed line, z = H). The f c mode is generated throughout the thickness 
of the disc at r = 1 but its energy rises from within the polytropic layer owing to wave channelling and becomes confined near the base 
of the atmosphere. The wave propagates without loss of flux tors 2.2 (Figure lOl dot-dashed line, left panel) before most of its energy 
is dissipated in two regions at the base of the atmosphere (white contour lines at r ?s 2.2). The dissipation in this example is partly 
numerical, as suggested by Figure [nj. If we were able to increase the resolution further, the wave would propagate a greater distance 
(c.f. different resolutions in the left panel of Figure |lo| ) and might agree even more closely with the linear prediction. 
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vertical energy distribution that is similar to the compress- 
ible case. Consequently, the process of energy concentration 
cannot be due to acoustic refraction. 

Figure ^l] shows an explicit comparison between the re- 
sults of the simulations and the predictions of acoustic re- 
fraction. In order to show both the wavefronts and the ver- 
tical extent of the wave, we plot a snapshot of the radial 
velocity, normalized by the RMS velocity, and multiplied by 
the mean wave energy density, from a simulation involving 
an f c mode in a polytropic disc. On top of the grey-scale, 
we trace the rays that would result from the refraction of an 
initially vertical wavefront at the resonant radius. The re- 
fracted wavefronts, which are orthogonal to the rays, would 
rapidly become severely tilted and the wave would be pre- 
dicted to propagate almost vertically into the isothermal 
atmosphere. However, the simulations show that the wave- 
fronts in fact remain nearly vertical as the wave propagates 
horizontally along the base of the disc atmosphere. The wave 
propagates without loss of flux to r ~ 2.2 before most of its 
energy is dissipated in two regions at the base of the at- 
mosphere (see Figure [To] for the radial energy flux). Figure 
pi] demonstrates that the linear theory of wave channelling 
accurately predicts the propagation of the wave until non- 
linear dissipation occurs. 

In summary, our non-linear hydrodynamic calculations 
clearly show that for axisymmetric waves the disc acts as 
a waveguide and that the behaviour of the waves is de- 
termined primarily by wave channelling and not refraction. 
Explanations based on simple refraction are incorrect. Ray 
tracing, if applied correctly, does offer a valid description of 
modes with a short vertical wavelength, i.e. those of large 
vertical mode number (see Figure 11 of Lubow & Pringle 
1993). However, such modes are unlikely to be excited by 
tidal resonant forcing and the f c mode, in particular, is ver- 
tically evanescent and has a vertical mode number of zero. 
As mentioned in Section 1, non-axisymmetric waves of low 
azimuthal wavenumber m ^ r/H are almost indistinguish- 
able from axisymmetric waves on the scale of a few H. We 
therefore expect that such waves will also propagate as if 
through a waveguide and be subject to similar wave chan- 
nelling in agreement with linear theory. Several papers (e.g. 
Lin et al. 1990a, b; Terquem 2001) have interpreted the be- 
haviour of low-m (e.g. m — 2) non-axisymmetric waves as 
being due to refraction. Given our results for m = waves, 
we conclude that the dominant behaviour of low-m non- 
axisymmetric waves is almost certainly determined by wave 
channelling and not refraction. 



The crest of the wave travels faster than the trough and the 
wave steepens as it propagates until shocks form. This is the 
principal effect leading to the dissipation of the f e mode in a 
vertically isothermal disc with 7 = 5/3 (Figure^). The situ- 
ation is not entirely straightforward because the wave Mach 
number is not uniform across the wavefront. In addition, 
the steepening may be accelerated by the gradual increase 
of wave amplitude as the wave propagates outwards into 
material of lower surface density while conserving its energy 
flux. On the other hand, the steepening may be temporarily 
postponed by the dispersive character of the wave close to 
the resonance. 

The r modes in the vertically isothermal disc are chan- 
nelled into an increasingly thin layer near the mid-plane and 
simultaneously develop very short radial wavelengths. These 
modes are therefore susceptible to viscous damping even in 
a nearly inviscid disc. In our simulations, these modes de- 
velop scales shorter than the grid spacing and are ultimately 
lost. 

In the polytropic disc, classical wave steepening does 
not occur in any simple sense because none of the waves be- 
haves like a plane acoustic wave. The f c mode, in particular, 
is highly dispersive, which tends to resist non- linear steep- 
ening. It undergoes rapid wave channelling to the base of 
the isothermal atmosphere and behaves like a surface grav- 
ity wave. The wave channelling enhances the amplitude of 
the wave by concentrating its energy into a smaller region 
of space. In earlier work (Lubow & Ogilvie 1998; Ogilvie & 
Lubow 1999) we showed that this effect is typically sufficient 
to amplify the wave to sonic velocities where steepening and 
dissipation are unavoidable. At the same time, the radial and 
vertical length-scales of the wave also decrease as it propa- 
gates outwards (Figure ^lj), making linear viscous damping 
a possible source of dissipation. Because wave channelling 
occurs over a distance that is almost independent of the 
thickness of the disc, the dissipation of waves with the same 
initial Mach number near the reasonance occurs at the same 
radius regardless of the disc's thickness. We have performed 
calculations with discs whose thickn esses vary by a factor of 
4 to demonstrate this effect (Section j.1.4 and Table [l]). The 
results of the simulations suggest that the f c mode in a poly- 
tropic disc dissipates energy at the base of the isothermal 
layer at a Mach number of approximately 0.2 (Figure |l2| ), 
which it acquires through the effect of wave channelling. The 
wave energy does not propagate vertically within the isother- 
mal atmosphere, contrary to the suggestion of Papaloizou & 
Lin (1995). 



8.3 Dissipation of waves 

In situations where waves in discs are excited by tidal forcing 
from an orbiting companion, the site at which a wave ulti- 
mately dissipates is of some importance, because the energy 
and angular momentum carried by the wave are transferred 
to the disc there. The simulations provide valuable informa- 
tion on the location and means of wave dissipation, matters 
about which we were previously able only to speculate (e.g. 
Lubow & Ogilvie 1998). 

In a nearly inviscid disc, dissipation can occur only if 
the wave develops very large velocity gradients. One way 
to achieve this is the classical non-linear steepening of an 
acoustic wave in a gas with 7 > 1 (e.g. Lighthill 1978). 



8.4 Outlook 

In future we hope to address some questions that remain 
unanswered in this paper. In particular, a more detailed 
study of the competition between linear dispersion and non- 
linear steepening of waves in discs would be valuable. Some 
important wave phenomena will require non-axisymmetric 
simulations. For example, a Keplerian disc supports slowly 
varying m = 1 density and bending waves that effect an ec- 
centric distortion and a warping of the disc, respectively. Un- 
like those studied in this paper, these waves can propagate 
over long distances without experiencing wave channelling. 
Finally, attention should be given to studying wave propa- 
gation in possibly more realistic models such as fully turbu- 
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lent discs and layered discs (Gammie 1996). These challenges 
provide ample opportunities for further developments. 



Yuan O, Cassen P., 1994, ApJ, 437, 338 



ACKNOWLEDGMENTS 

We are grateful to Jim Stone for providing the ZEUS-2D 
code and for many useful discussions. We also gratefully 
acknowledge support from NASA grants NAG5-4310 and 
NAG5-10732, the STScI visitor program, and the Institute 
of Astronomy visitor programme. Some of the calculations 
discussed in this paper were performed on the GRAND com- 
puter. GIO acknowledges the support of the Royal Society 
through a University Research Fellowship. 



REFERENCES 

Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical 

Functions, Dover, New York 
Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651 
Bell K. R., Cassen P. M., Klahr H. H., Henning Th., 1997, ApJ, 

486, 372 

Cassen P., Woolum D. S. 1996, ApJ, 472, 789 

Gammie C. F., 1996, ApJ, 457, 355 

Goldreich P., Tremaine S. 1979, ApJ, 233, 857 

Goldreich P., Tremaine S. 1980, ApJ, 241, 425 

Jackson J. D., 1962, Classical Electrodynamics, Wiley, New York 

Kato S., Fukue J., 1980, PASJ, 32, 377 

Kato S., 2001, PASJ, 53, 1 

Korycansky D. C, Pringle J. E., 1995, MNRAS, 272, 618 
Landau L. D., Lifshitz E. M., 1960, Electrodynamics of Continu- 
ous Media, Pergamon, New York 
Larson R. G., 1990, MNRAS, 243, 588 

Lighthill M. J., 1978, Waves in Fluids, Cambridge Univ. Press, 
Cambridge 

Lin D. N. C, Papaloizou J. C. B., 1979, MNRAS, 186, 799 

Lin D. N. C, Papaloizou J. C. B., 1993, in Levy E. H., Lunine 

J. I., eds, Protostars and Planets III, Univ. Arizona Press, 

Tucson, p. 749 

Lin D. N. C, Papaloizou J. C. B., Savonije G. J., 1990a, ApJ, 

364, 326 

Lin D. N. C, Papaloizou J. C. B., Savonije G. J., 1990b, ApJ, 

365, 748 

Loska Z., 1986, Acta Astron., 36, 43 
Lubow S. H., 1981, ApJ, 245, 274 
Lubow S. H., 1991, ApJ, 381, 259 

Lubow S. H., Artymowicz P., 1996, in Wijers R. A. M. J., Davies 
M. B., Tout C. A., eds, Evolutionary Process in Binary Stars, 
Kluwer, Dordrecht, p. 53 
Lubow S. H., Ogilvie G. I., 1998, ApJ, 504, 983 
Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360 
Nowak M. A., Wagoner, R. V., 1991, ApJ, 378, 656 
Ogilvie G. I., 1998, MNRAS, 297, 291 
Ogilvie G. I., 2002, MNRAS, in press 
Ogilvie G. I., Lubow S. H., 1999, ApJ, 515, 767 
Okazaki A. T., Kato S., Fukue J., 1987, PASJ, 39, 457 
Papaloizou J. C. B., Lin D. N. C, 1995, ARA&A, 33, 505 
Papaloizou J. C. B., Nelson R. P., Masset F., 2001, A&A, 366, 
263 

Shu F. H., Cuzzi J. N., Lissauer J. J., 1983, Icarus, 53, 185 
Shu F. H., Yuan C, Lissauer J. J., 1985, ApJ, 291, 356 
Stone J. M., Norman M. L., 1992, ApJS, 80, 753 
Terquem C. E. J. M. L. J., 2001, in Mathieu R. D., Zinnecker H., 

eds, The Formation of Binary Stars, ASP Conf. Ser., vol. ???, 

p.??? 

von Neumann J., Richtmyer R. D., 1950, J. Appl. Phys., 21, 232 



APPENDIX A: DISC MODELS 

Al A polytropic disc with an isothermal 
atmosphere 

We consider a differentially rotating fluid with density 
p(R,z), pressure p(R,z) and angular velocity £l(R, z), re- 
ferred to cylindrical polar coordinates (R,<f),z). The equa- 
tions of equilibrium are 



ROT = 
= 



1 dp <9<E> 
pdR ~ OR' 
1 dp d& 

p dz dz' 



where the gravitational potential is 
$(R,z) = -GM(R 2 + z 2 y 1/2 . 



(Al) 
(A2) 

(A3) 



Let the pressure and density be related by (cf. Lin, Pa- 
paloizou & Savonije 1990) 



2 
C p 



l/n 
~ 1 + J 

P> 



(A4) 



where n is a positive constant, while c(R) and p s (R) are pos- 
itive functions. At high densities p 3> p s , the gas becomes 
vertically polytropic, with polytropic index n. At low den- 
sities p <C p s , the gas becomes vertically isothermal, with 
sound speed c. The transitional density p s defines the ap- 
proximate surface of the disc. 
Define the pseudo-enthalpy 



h — c 



l/n 



■ In I — 1 - 1 



Then it follows that 
dp 
P 

where 



= dh + FdR, 
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(A7) 



Consider first the vertical equilibrium at each R sepa- 
rately. Equation (A2) becomes 







dh _ 9£ 
dz dz ' 
with solution 



(A8) 



ft =-* + /(«). (A9) 
By considering the mid-plane z — we determine 

l/n 



+ln(^|-l 



GM 



(A10) 



Waves in accretion discs 25 



where pc(R) is the mid-plane density. In general, therefore, 

l/n / \ l/n 



P_\ _ I Pc 

Ps J \ps 

GM 



+ ln 



GM 



(All) 



Given the functions p c (R), p s (R) and c(R), this equation can 
be solved numerically to dete rmine p at any given point, and 
p follows from equation (A4). 



Now consider the radial equilibrium. Equation (Al) be- 
comes 



dh 
OR 



■F + 



dR' 



Thus, from equation (A9), 
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(A13) 



which determines the angular velocity at every point. 

Let us assume that the disc is thin and that p c ^> p s so 
that there is a well-defined polytropic layer with an isother- 
mal atmosphere outside. The scale-height of the isothermal 
layer is 

while the thickness of the polytropic layer is 
l/2n 



,1/2 



H. 



(A15) 



If we take c oc R~ 1,>2 and p c /p s = constant then Hi oc 
H p oc R so that the scale-height of the disc increases linearly 
with radius. We define f3 = p s /pc and e to be the nominal 
value of Hp/R. Thus 



-l/2fll/2n 



/GM \ 1/2 



V R J 



(A16) 



c = 2 
and 

p s = pp c . (A17) 

Finally, we take the profile of p c (R) to be a power law but 
tapered to give inner and outer edges, 

' R — R\ \ 1 . / R2 — R 



Pc = R~ 



c 1 , /J 

+ — tann I - 



1 

+ — tanh 



(A18) 



where S <C 1 is the factor by which the density is reduced 
outside the interval i?i <i?<i?2 occupied by the disc, and 
wi and W2 are the widths of the inner and outer tapers. The 
parameters of the model are then n,e,/3, a, 5, Ri , R2 , Wi , W2 . 

For our standard polytropic discs, we choose G = 1, 
M = 1, n = 1.5, (5 = 0.01, a = 1.5, S = 0.01, R 2 = 3.0, wi = 



Rit, and W2 = -R2C The value of Ri = 0.2 for all but the 
highest resolution calculations, for which _Ri = 0.6. Results 
are given for three different values of the disc thickness, e = 
Hp/R = 0.05, 0.1, 0.2. 

A2 A vertically isothermal disc 

In the case of a vertically isothermal disc, we take 
V = c 2 p, (A19) 
P X 



Pc 



(A20) 



with c = c(R), p = Pc(R)- Then the analysis proceeds as 
before, but now with 
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The angular velocity is given by 
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To obtain a scale-height that increases linearly with radius 
we take 



(A26) 



where, as before, e is a small constant, so that the isothermal 
scale-height H — c/fi satisfies H/R ~ e. As in the p olytropic 
case, the profile of p c (R) is given by equation |A1S| . 

We use the same values of G, M, a, 5 for the vertically 
isothermal discs as we use for the polytropic discs, but take 
Ri — 0.5, R2 = 10.0, «ii = Rie and W2 = ifee, while n and 
are not used. 



APPENDIX B: RESONANT LAUNCHING OF 
AXISYMMETRIC WAVES IN LINEAR THEORY 

Bl The f c mode 

The launching of non-axisymmetric waves at a Lindblad 
resonance in a three-dimensional disc has been analysed 
by Lubow & Ogilvie (1998). The total torque exerted at 
the resonance is identical to that determined by Goldre- 
ich & Tremaine (1979) using a two-dimensional model that 
neglects the vertical structure of the disc. However, in a 
three-dimensional disc the torque is partitioned into differ- 
ent wave modes. In an inviscid disc, the total torque is equal 
to the sum of the radial fluxes of angular momentum of the 
launched waves, averaged over time, integrated azimuthally 
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and vertically, and reckoned in the direction of propagation 
at some distance from the resonance. 

In the case of axisymmetric forcing no torque is exerted. 
However, the work done at the resonance can still be deter- 
mined from the formula of Goldreich & Tremaine (1979) in 
an appropriate limit, using the fact that the energy flux of 
each wave is equal to the angular momentum flux multiplied 
by the angular pattern frequency. For a Keplerian disc the 
rate at which work is done is equal to 



T ■ = 



n 2 R 2 A 2 T, 
3uj 



(Bl) 



where E = J p dz is the surface density, and the applied 
radial acceleration is cir — A(R) exp(— iuit). All quantities 
are evaluated at the location of the resonance. The Lindblad 
resonance condition ui = k corresponds to u — Q in a Kep- 
lerian disc. In an inviscid disc, T is equal to the total energy 
flux in the launched waves, i.e. the sum of their radial en- 
ergy fluxes, averaged over time, integrated azimuthally and 
vertically, and reckoned in the direction of propagation at 
some distance from the resonance. 

The fraction of the energy flux imparted to any wave 
mode is given by equation (45) of Lubow & Ogilvie (1998), 



/ 



pu'ji dz 




P\u R \ 



' dz, 



(B2) 



where u' R (z) is the eigenfunction of the radial velocity per- 
turbation for the mode considered, and the integrals are over 
the full vertical extent of the disc. Again, all these quanti- 
ties are evaluated at the resonance. The imparted fraction 
therefore depends on an overlap integral between the applied 
forcing and the modes of the disc. 

Lubow & Ogilvie (1998) found that, in a polytropic disc, 
more than 95% of the flux is taken up by the f° mode. The 
remainder goes mainly into the r° mode. We have repeated 
the calculation for the standard polytropic model adopted 
in the present paper, including the effects of the isothermal 
atmosphere and the vertical boundary conditions imposed 
in the simulations. We find that 98.3% and 1.6% should go 
into the f e and r 1 ; modes, respectively. We recall that the 
r modes propagate away from the resonance in the opposite 
direction to the f c mode. The vertical boundaries are suffi- 
ciently distant that they have no effect on the modes or on 
the flux fractions. 

In a vertically isothermal disc the eigenfunctions are 
known analytically from the work of Lubow & Pringle 
(1993). The radial velocity perturbation for the two- 
dimensional (or f c ) mode is of the form^j 



Un oc exp 



7 



z 

2H 2 



(B3) 



while, of course, p oc exp(— z 2 /2H 2 ). We then find for the 
fraction imparted to the two-dimensional (or f c ) mode 



/2D = [ 7 (2 - 7)] 



1/2 



(B4) 



Only in the case of isothermal perturbations (7 = 1) does the 
forcing match perfectly to the two-dimensional mode. When 

1 There is a sign error in the argument of the exponential in the 
first and second lines of equation (37) of Lubow & Pringle (1993). 



7 = 5/3, as in the present paper, /2D = v5/3 « 74.5%. A 
further 6.0% goes into the rf mode, and yet smaller frac- 
tions into higher-order r modes. Interestingly, 18.9% of the 
expected flux is not accounted for by previously documented 
modes, which evidently do not form a complete set. 

When we calculate the eigenfunctions of the modes at 
a Lindblad resonance in a vertically isothermal disc using 
the vertical boundary conditions imposed in the simula- 
tions, the flux fractions imparted to the two-dimensional 
and r modes agree with the numbers given above, indicat- 
ing that the boundaries are sufficiently distant not to affect 
the modes. At the same time, we discover a further sequence 
of modes that fully account for the 'missing' 18.9% of the 
expected flux. These modes propagate on the same side of 
the resonance as the f c mode but are strongly affected by the 
boundaries. They become spectrally dense in the limit that 
the boundaries are removed to a great height, and disappear 
in the limit 7 — > 1. Most of the wave energy of these modes 
is contained away from the mid-plane, near the boundaries. 

We have identified these additional modes as g modes 
that are artificially confined by the boundaries. In a verti- 
cally isothermal disc with 7 > 1, the buoyancy frequency 
is proportional to z and therefore increases indefinitely with 
height. In the absence of vertical boundaries, internal gravity 
waves launched at a Lindblad resonance may propagate with 
a positive but ever diminishing vertical group velocity and 
never reach a turning point. As a result, they can never sat- 
isfy a standing wave condition and form a vertically confined 
g mode. However, they constitute a continuous spectrum of 
outgoing waves and, as such, receive a non-zero share of the 
total energy flux. 

High in the atmosphere of a real disc, the adiabatic 
exponent 7 = 1, since the compressional energy of the wave 
is radiated away in less than a wave period. In that case, 
we expect that some g modes will be excited and will be 
vertically confined. The individual g modes that are excited 
should take up the required energy flux in a manner that is 
independent of the location of any vertical boundaries. 

If we start from our standard polytropic model and con- 
tinuously reduce the effective optical thickness towards zero 
by increasing the parameter j3 so that the isothermal atmo- 
sphere acquires more mass, the flux fraction for the f e mode 
decreases continuously from 98.3% towards 74.5%. 



B2 The p\ mode 

The pi mode is launched at a vertical resonance of the type 
analysed by Lubow (1981). In the case of axisymmetric forc- 
ing the resonance condition is uj = (7 + l) 1//2 n z , where Q, z is 
the vertical oscillation frequency, equal to Q is a Keplerian 
disc. Near the resonance, the pf mode involves a predomi- 
nantly vertical velocity perturbation u' z oc z, corresponding 
to a 'breathing' motion of the disc, and is excited efficiently 
by an applied vertical acceleration a z = A(R)z exp(— iut). 
We give below an informal derivation of the wave equation 
satisfied by the p\ mode near the resonance, in order to 
predict the energy flux in the launched wave. A rigorous 
derivation can be given, if desired, using asymptotic expan- 
sions. 

The exact linearized equations governing axisymmetric 
Eulerian perturbations of the disc in the presence of vertical 
forcing of the above type are 
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In a thin, Keplerian disc we may take Q, = (GM/R 3 ) 1 ' 2 , 
neglecting its vertical variation, and set dp/dz = —pQ 2 z for 
vertical hydrostatic equilibrium. For a disturbance with ra- 
dial wavelength short compared to R, we may also ne glec t 
the terms involving radial derivatives of p in equation (B5), 
of p and of R in equation (B£), and of p and of R in equa- 
tion (B9). With these approximations we obtain, after an 



integration by parts, 
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The quantity in square brackets vanishes at the res- 
onance. If X = R — R ICS is the radial distance from the 
resonance, we may apply the linear approximation 



o; 2 -( 7 + l)fi 2 



n 2 

3(7+l)^X 



(Bll) 



The leading approximation to the p° mode near the reso- 
nance is 

u' z = iiozY(X), 



= (>+'l)™. 



(B12) 



where Y(X) is a function to be determined. The horizontal 
components of the equation of motion then determine 

dp} dY_ 
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The integral relation (BIO) then yields the ordinary differ- 
ential equation 
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where the coefficients 
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are all to be evaluated at the location of the resonance. By 
means of the change of variables 



Y(X) 
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equation (B14) is transformed into the inhomogeneous Airy 
equation, 

as also appears in the analysis of Lindblad resonances. The 
desired solution, representing a wave emitted from the res- 
onance and travelling into X > with no incoming compo- 
nent, is 



y = -7r [Gi(-x) + i Ai(-a 



(B18) 



where Ai and Gi are homogeneous and inhomogeneous 
Airy functions (e.g. Abramowitz & Stegun 1965). From the 
asymptotic form for large positive x, 
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we determine the radial energy flux (averaged in time, and 
integrated azimuthally and vertically) in the wave at some 
distance from the resonance, 



nR J u' R *p dz, 
n 2 R 2 A 2 
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This is directly analogous to equation (Bl) for a Lindblad 
resonance. 



B3 The f°/rf mode 

Generically, the f° mode, which effects a corrugation or tilt 
of the disc, is launched at a vertical resonance. In the case of 
axisymmetric forcing the resonance condition is uj = £l z . In 
a Keplerian disc this coincides with the Lindblad resonance 
and the f°/r° mode is launched at a hybrid vertical/Lindblad 
resonance that has not been discussed in the literature. 

Near the resonance, the f°/r° mode involves a vertical 
velocity perturbation u' z independent of z and horizontal 
velocity perturbations u' R , ti'^ a z of comparable magnitude. 
It is excited efficiently by an applied vertical acceleration 
a z = A(R) exp(— iujt). We give below an informal derivation 
of the wave equation satisfied by the f°/r° mode near the 
resonance. Again, a rigorous derivation can be given using 
asymptotic expansions. 

The equations governing the perturbations are the same 
as in the previous subsection except that the vertical accel- 
eration no longer contains the factor z. We now obtain the 
integral relation 

iuA J p dz = {to 2 - n 2 ) J pu' z dz- j |2 dz, (B21) 

in which the quantity in brackets vanishes at the resonance. 

If X = R — R ICS is the radial distance from the reso- 
nance, we may apply the linear approximation 
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The leading approximation to the f°/r° mode near the res- 
onance is 

u' z = iuY(X), 
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The horizontal components of the equation of motion then 
determine 
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The denominator of this expression may be r epresented by 
the same linea r approximation, equation (B22). The integral 



relation (B21) then yields the ordinary differential equation 
d 



where the coefficients 
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are all to be evaluated at the location of the resonance. By 
means of the change of variables 
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equation 0B25I) is transformed into 
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Unlike equation (B17), this equation allows wave propaga- 
tion on both sides of the resonance, x > and x < 0. The 



complementary functions of equation (B28) are cos(^x 



and sin(^a; 2 ), and the general solution of the inhomogeneous 
equation can therefore be expressed in terms of Fresnel inte- 
grals. The desired solution, representing one wave travelling 
into X > and another wave travelling into X < 0, with no 
incoming component, is 
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From the asymptotic form for large positive x, 
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we determine the energy flux in the wave at some distance 
from the resonance on the positive side, 



T — nR J u' R *p'dz, 
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This is exactly half the value given by equation (Bl) for a 



Lindblad resonance. From the symmetry property y(—x) = 
—y*(x), it follows that the wave emitted into x < carries 
an equal and oppositely directed energy flux. Therefore the 
total energy flux is identical to that for a Lindblad reso- 
nance, although of course A is the amplitude of the vertical, 
not horizontal, acceleration. 
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